Optimized functions to compute projection of a point on a line?

Question

Consider the following problem:

My question is the following: how to optimize the following independent functions:

// Computation of the coordinates of P
inline std::array<double, 3> P(const std::array<double, 3>& A, 
                               const std::array<double, 3>& B,
                               const std::array<double, 3>& M)
{
    // The most inefficient version in the world (to be verified)
    std::array<double, 3> AB = {B[0]-A[0], B[1]-A[1], B[2]-A[2]};
    std::array<double, 3> AM = {M[0]-A[0], M[1]-A[1], M[2]-A[2]};
    double norm = std::sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]);
    double dot = AB[0]*AM[0]+AB[1]*AM[1]+AB[2]*AM[2];
    double d1 = dot/norm;
    std::array<double, 3> AP = {AB[0]/d1, AB[1]/d1, AB[2]/d1};
    std::array<double, 3> P = {AP[0]-A[0], AP[1]-A[1], AP[2]-A[2]};
    return P;
}

// Computation of the distance d0
inline double d0(const std::array<double, 3>& A, 
                 const std::array<double, 3>& B,
                 const std::array<double, 3>& M)
{
    // The most inefficient version in the world (to be verified)
    std::array<double, 3> AB = {B[0]-A[0], B[1]-A[1], B[2]-A[2]};
    std::array<double, 3> AM = {M[0]-A[0], M[1]-A[1], M[2]-A[2]};
    double norm = std::sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]);
    double dot = AB[0]*AM[0]+AB[1]*AM[1]+AB[2]*AM[2];
    double d1 = dot/norm;
    std::array<double, 3> AP = {AB[0]/d1, AB[1]/d1, AB[2]/d1};
    std::array<double, 3> P = {AP[0]-A[0], AP[1]-A[1], AP[2]-A[2]};
    std::array<double, 3> MP = {P[0]-M[0], P[1]-M[1], P[2]-M[2]};
    double d0 = std::sqrt(MP[0]*MP[0]+MP[1]*MP[1]+MP[2]*MP[2]);
    return d0;
}

// Computation of the distance d1
inline double d1(const std::array<double, 3>& A, 
                 const std::array<double, 3>& B,
                 const std::array<double, 3>& M)
{
    // The most inefficient version in the world (to be verified)
    std::array<double, 3> AB = {B[0]-A[0], B[1]-A[1], B[2]-A[2]};
    std::array<double, 3> AM = {M[0]-A[0], M[1]-A[1], M[2]-A[2]};
    double norm = std::sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]);
    double dot = AB[0]*AM[0]+AB[1]*AM[1]+AB[2]*AM[2];
    double d1 = dot/norm;
}

// Computation of the distance d2
inline double d2(const std::array<double, 3>& A, 
                 const std::array<double, 3>& B,
                 const std::array<double, 3>& M)
{
    // The most inefficient version in the world (to be verified)
    std::array<double, 3> AB = {B[0]-A[0], B[1]-A[1], B[2]-A[2]};
    std::array<double, 3> AM = {M[0]-A[0], M[1]-A[1], M[2]-A[2]};
    double norm = std::sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]);
    double dot = AB[0]*AM[0]+AB[1]*AM[1]+AB[2]*AM[2];
    double d1 = dot/norm;
    double d2 = norm-d1;
    return d2;
}

So that each function will be as much optimized as possible ? (I will execute these functions billion times).

Answer 1

From algorithm point of view, you can calculate projection of vector to another vector not using SQRT call

the pseudocode from here http://www.euclideanspace.com/maths/geometry/elements/line/projections/

// projection of vector v1 onto v2
inline vector3 projection( const vector3& v1, const vector3& v2 ) {
    float v2_ls = v2.len_squared();
        return v2 * ( dot( v2, v1 )/v2_ls );

}

where dot() is a dot product of two vectors and len_squared is the dot product of vector with self.

NOTE: Try to pre calculate inverse of v2_ls before main loop, if possible.