Taylor McLaughlin Series to estimate the distance of two points

Question

Distance from point to point: dist = sqrt(dx * dx + dy * dy); But sqrt is too slow and I can't accept that. I found a method called Taylor McLaughlin Series to estimate the distance of two points on the book. But I can't comprehend the following code. Thanks for anyone who helps me.

#define MIN(a, b) ((a < b) ? a : b)
int FastDistance2D(int x, int y)
{
    // This function computes the distance from 0,0 to x,y with 3.5% error
    // First compute the absolute value of x, y
    x = abs(x);
    y = abs(y);

    // Compute the minimum of x, y
    int mn = MIN(x, y);

    // Return the distance
    return x + y - (mn >> 1) - (mn >> 2) + (mn >> 4);
}

I have consulted related data about McLaughlin Series, but I still can't comprehend how the return value use McLaughlin Series to estimate the value. Thanks for everyone~

Answer 1

This task is almost duplicate of another one: Very fast 3D distance check?

And there was link to great article: http://www.azillionmonkeys.com/qed/sqroot.html

In the article you can find different aproaches for approximation of root. For example maybe this one is suitable for you:

int isqrt (long r) {
    float tempf, x, y, rr;
    int is;

    rr = (long) r;
    y = rr*0.5;
    *(unsigned long *) &tempf = (0xbe6f0000 - *(unsigned long *) &rr) >> 1;
    x = tempf;
    x = (1.5*x) - (x*x)*(x*y);
    if (r > 101123) x = (1.5*x) - (x*x)*(x*y);
    is = (int) (x*rr + 0.5);
    return is + ((signed int) (r - is*is)) >> 31;
}

If you can calculate root operation fast, then you can calculate distance in regular way:

return isqrt(a*a+b*b)

And one more link: http://www.flipcode.com/archives/Fast_Approximate_Distance_Functions.shtml

u32 approx_distance( s32 dx, s32 dy )
{
   u32 min, max;

   if ( dx < 0 ) dx = -dx;
   if ( dy < 0 ) dy = -dy;

   if ( dx < dy )
   {
      min = dx;
      max = dy;
   } else {
      min = dy;
      max = dx;
   }

   // coefficients equivalent to ( 123/128 * max ) and ( 51/128 * min )
   return ((( max << 8 ) + ( max << 3 ) - ( max << 4 ) - ( max << 1 ) +
            ( min << 7 ) - ( min << 5 ) + ( min << 3 ) - ( min << 1 )) >> 8 );
}