Jump to content
  • Advertisement
Sign in to follow this  
Jorn Bloodyscalpel

Could you tell me how the code works ?

This topic is 4654 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
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!