It is probably not the quickest reply, but it surely’s positively one of many cleanest:
You can convert to homogenous coordinates in an effort to use Cramer’s rule cleanly for every barycentric coordinate.
Please ignore the unusual syntax, it’s a customized preprocessor I wrote over OpenCL C.

VEC3 barycentric(VEC3 P, Triangle3f tri){

    VEC3 A = VEC3::diff(_VEC3(tri.A),P);
    VEC3 B = VEC3::diff(_VEC3(tri.B),P);
    VEC3 C = VEC3::diff(_VEC3(tri.C),P);
    
    float Ax = A.unwrap()..x;
    float Ay = A.unwrap()..y;
    float Az = A.unwrap()..z;
    
    float Bx = B.unwrap()..x;
    float By = B.unwrap()..y;
    float Bz = B.unwrap()..z;
    
    float Cx = C.unwrap()..x;
    float Cy = C.unwrap()..y;
    float Cz = C.unwrap()..z;
    
    float Px = 0.0;
    float Py = 0.0;
    float Pz = 0.0;
    
    MAT4 M = __MAT4(
    
        Ax,Bx,Cx,1,
        Ay,By,Cy,1,
        Az,Bz,Cz,1,
        1, 1, 1, 0
    
    );
    
    float d = MAT4::det(M).scalarUnwrap();
    
    if(fabs(d)<DET_EPSILON){
        return VEC3::zero();
    }
    
    MAT4 M1 = __MAT4(
    
        Px,Bx,Cx,1,
        Py,By,Cy,1,
        Pz,Bz,Cz,1,
        1, 1, 1, 0
    
    );
    
    MAT4 M2 = __MAT4(
    
        Ax,Px,Cx,1,
        Ay,Py,Cy,1,
        Az,Pz,Cz,1,
        1, 1, 1, 0
    
    );
    
    MAT4 M3 = __MAT4(
    
        Ax,Bx,Px,1,
        Ay,By,Py,1,
        Az,Bz,Pz,1,
        1, 1, 1, 0
    
    );
    
    return VEC3::fromcoords(
        MAT4::det(M1).scalarUnwrap() / d,
        MAT4::det(M2).scalarUnwrap() / d,
        MAT4::det(M3).scalarUnwrap() / d
    );

}