• Advertisement
Sign in to follow this  

Could you tell me how the code works ?

This topic is 4412 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

I have got some code which calculates the length of a vector, but i don't understand it. The code is showed below: ======================================================== int Fast_Distance_2D(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)); } // end Fast_Distance_2D /////////////////////////////////////////////////////////////////////////////// float Fast_Distance_3D(float fx, float fy, float fz) { // this function computes the distance from the origin to x,y,z int temp; // used for swaping int x,y,z; // used for algorithm // make sure values are all positive x = fabs(fx) * 1024; y = fabs(fy) * 1024; z = fabs(fz) * 1024; // sort values if (y < x) SWAP(x,y,temp) if (z < y) SWAP(y,z,temp) if (y < x) SWAP(x,y,temp) int dist = (z + 11 * (y >> 5) + (x >> 2) ); // compute distance with 8% error return((float)(dist >> 10)); } // end Fast_Distance_3D ====================================================== The author said it's base on Taylor Series. I know Taylor Series, but I don't know how the code works. Could you help me? Thank you! I'm not a native English speaker, and I hope my words are clear enough.

Share this post


Link to post
Share on other sites
Advertisement
It is simply an approximation. In the 2D version, the function d = sqrt(a2 + b2) is approximated by the function d = a + .375b(where a is the length of the longest side and b is the length of the shortest side). Using shortest and longest reduces the error. The 3D version is similar.

Here is a graph of the two functions (a is constant):

Image Hosted by ImageShack.us


[Edited by - JohnBolton on January 20, 2006 3:49:09 PM]

Share this post


Link to post
Share on other sites
With all those branches that 'fast approximation' is probably slower than doing it the normal way (and without the 8% error). Maybe someone less lazy than me can whip up a quick benchmark.

Share this post


Link to post
Share on other sites
i think it's called the manhattan distance.

it uses maclaurin/taylor series to approximate, however i don't know how it's derived since it has two variables and i'm very curious. if anyone can explain?

Share this post


Link to post
Share on other sites
Quote:
Original post by OrangyTang
With all those branches that 'fast approximation' is probably slower than doing it the normal way (and without the 8% error). Maybe someone less lazy than me can whip up a quick benchmark.


I did a quick test, I got the following results:
Fast_Distance_3D: 10.3533 secs, 2^32/10 iterations
sqrt(x*x+y*y+z*z): 7.68405 secs, 2^32/10 iterations
Fast_Distance_2D: 40.7154 secs, 2^32 iterations
int_sqrt(x*x+y*y): 1.39683E-006 secs, 2^32 iterations

For integer sqrt I used the following function:

unsigned int isqrt(unsigned int input)
{
unsigned int rest = input;
unsigned int step = (1<<30);
unsigned int result = 0;

if (step <= rest)
result |= (1<<15), rest -= step;
step = (result << 15) + (1<<28);
if (step <= rest)
result |= (1<<14), rest -= step;
step = (result << 14) + (1<<26);
if (step <= rest)
result |= (1<<13), rest -= step;
step = (result << 13) + (1<<24);
if (step <= rest)
result |= (1<<12), rest -= step;
step = (result << 12) + (1<<22);
if (step <= rest)
result |= (1<<11), rest -= step;
step = (result << 11) + (1<<20);
if (step <= rest)
result |= (1<<10), rest -= step;
step = (result << 10) + (1<<18);
if (step <= rest)
result |= (1<<9), rest -= step;
step = (result << 9) + (1<<16);
if (step <= rest)
result |= (1<<8), rest -= step;
step = (result << 8) + (1<<14);
if (step <= rest)
result |= (1<<7), rest -= step;
step = (result << 7) + (1<<12);
if (step <= rest)
result |= (1<<6), rest -= step;
step = (result << 6) + (1<<10);
if (step <= rest)
result |= (1<<5), rest -= step;
step = (result << 5) + (1<<8);
if (step <= rest)
result |= (1<<4), rest -= step;
step = (result << 4) + (1<<6);
if (step <= rest)
result |= (1<<3), rest -= step;
step = (result << 3) + (1<<4);
if (step <= rest)
result |= (1<<2), rest -= step;
step = (result << 2) + (1<<2);
if (step <= rest)
result |= (1<<1), rest -= step;
step = (result << 1) + (1<<0);
if (step <= rest)
result++;

if (input > result*result+result) result++; //for correct rounding
return result;
}




I did the test on an AMD Athlon64 3200+ with 1024 MB DDR RAM in Release mode with MSVC2005. Also I used QueryPerformanceCounter to get time.

Here is the complete test code:

#include <iostream>
#include <cmath>
#include <limits>
#include <windows.h>

#define SWAP(i,j,k) k = i; i = j; j = k
#define MIN(x,y) x>y?x:y

unsigned int isqrt(unsigned int input)
{
unsigned int rest = input;
unsigned int step = (1<<30);
unsigned int result = 0;

if (step <= rest)
result |= (1<<15), rest -= step;
step = (result << 15) + (1<<28);
if (step <= rest)
result |= (1<<14), rest -= step;
step = (result << 14) + (1<<26);
if (step <= rest)
result |= (1<<13), rest -= step;
step = (result << 13) + (1<<24);
if (step <= rest)
result |= (1<<12), rest -= step;
step = (result << 12) + (1<<22);
if (step <= rest)
result |= (1<<11), rest -= step;
step = (result << 11) + (1<<20);
if (step <= rest)
result |= (1<<10), rest -= step;
step = (result << 10) + (1<<18);
if (step <= rest)
result |= (1<<9), rest -= step;
step = (result << 9) + (1<<16);
if (step <= rest)
result |= (1<<8), rest -= step;
step = (result << 8) + (1<<14);
if (step <= rest)
result |= (1<<7), rest -= step;
step = (result << 7) + (1<<12);
if (step <= rest)
result |= (1<<6), rest -= step;
step = (result << 6) + (1<<10);
if (step <= rest)
result |= (1<<5), rest -= step;
step = (result << 5) + (1<<8);
if (step <= rest)
result |= (1<<4), rest -= step;
step = (result << 4) + (1<<6);
if (step <= rest)
result |= (1<<3), rest -= step;
step = (result << 3) + (1<<4);
if (step <= rest)
result |= (1<<2), rest -= step;
step = (result << 2) + (1<<2);
if (step <= rest)
result |= (1<<1), rest -= step;
step = (result << 1) + (1<<0);
if (step <= rest)
result++;

if (input > result*result+result) result++; //for correct rounding
return result;
}

int Fast_Distance_2D(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));

} // end Fast_Distance_2D

float Fast_Distance_3D(float fx, float fy, float fz)
{
int temp;
int x,y,z; // used for algorithm

// make sure values are all positive
x = static_cast<int>(fabs(fx) * 1024.0F);
y = static_cast<int>(fabs(fy) * 1024.0F);
z = static_cast<int>(fabs(fz) * 1024.0F);

// sort values
if (y < x) SWAP(x,y,temp);

if (z < y) SWAP(y,z,temp);

if (y < x) SWAP(x,y,temp);

int dist = (z + 11 * (y >> 5) + (x >> 2) );

// compute distance with 8% error
return((float)(dist >> 10));
}
float Distance_3D(float fx, float fy, float fz)
{
return sqrtf(fx*fx+fy*fy+fz*fz);
}
int Distance_2D(int fx, int fy)
{
return isqrt(fx*fx+fy*fy);
}

#undef max

int main()
{
const unsigned int num_iter = std::numeric_limits<unsigned int>::max()/10;
const unsigned int num_iter2 = std::numeric_limits<unsigned int>::max();

float result = 0.0F;
int result2 = 0;

__int64 freq;
__int64 time0;
__int64 time1;
QueryPerformanceFrequency( reinterpret_cast<LARGE_INTEGER*>(&freq));
QueryPerformanceCounter( reinterpret_cast<LARGE_INTEGER*>(&time0));

for(unsigned int i = 0;
i < num_iter;
++i)
{
result += Fast_Distance_3D( result, result, result);
}

QueryPerformanceCounter( reinterpret_cast<LARGE_INTEGER*>(&time1));

std::cout << "Time to use \"fast\" distance function: " <<
static_cast<double>(time1 - time0)/static_cast<double>(freq) << std::endl;

QueryPerformanceCounter( reinterpret_cast<LARGE_INTEGER*>(&time0));

for(unsigned int i = 0;
i < num_iter;
++i)
{
result += Distance_3D( result, result, result);
}

QueryPerformanceCounter( reinterpret_cast<LARGE_INTEGER*>(&time1));

std::cout << "Time to use distance function: " <<
static_cast<double>(time1 - time0)/static_cast<double>(freq) << std::endl;


QueryPerformanceCounter( reinterpret_cast<LARGE_INTEGER*>(&time0));

for(unsigned int i = 0;
i < num_iter2;
++i)
{
result2 += Fast_Distance_2D( result2, result2);
}

QueryPerformanceCounter( reinterpret_cast<LARGE_INTEGER*>(&time1));

std::cout << "Time to use \"fast\" 2d distance function: " <<
static_cast<double>(time1 - time0)/static_cast<double>(freq) << std::endl;

QueryPerformanceCounter( reinterpret_cast<LARGE_INTEGER*>(&time0));

for(unsigned int i = 0;
i < num_iter2;
++i)
{
result2 += Distance_2D( result2, result2);
}

QueryPerformanceCounter( reinterpret_cast<LARGE_INTEGER*>(&time1));

std::cout << "Time to use 2d distance function: " <<
static_cast<double>(time1 - time0)/static_cast<double>(freq) << std::endl;

std::cout << "Garbage:" << result <<std::endl;

std::cin.get();
return 0;
}


Share this post


Link to post
Share on other sites
Well, I know, the discussion about sqrt has nothing new (I have checked all the threads before). I just feel curious in how the formula was derived (though i prefer carmack sqrt in fact). I think the method used in the code i posted may get well used in other calculating.

To CTar:

The isqrt function you posted i have found in glu lib before, but i don't understand it, too. Could you or anyone else do me a favor to explain it ?

Thank you very much !

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement