Calculating square roots.

[hamster_nz]

Here is an interesting way to calculate the square root of an integer that would be suited to FPGA implementation. It only uses bit shifts and addition/subtraction so might be of use to micro users

Code: [Select]

/*****************************************
* my_sqrt.c: Implmenting an integer square
*            root using only bit-shifts
*
* Author: Mike Field <hamster@snap.net.nz>
*
*****************************************/
#include <stdio.h>
#include <math.h>

unsigned my_sqrt(unsigned a) {
  unsigned s = 0; /* Where the result is assembled */
  int      q = 15;  /* Bit we are working on */

  while(q >= 0) {
    /* 'change' represents how much setting bit
     * 'q' in 's' will change the value of s^2
     *
     * It is (2^q)*(2^q) + s*(2^q) + s*(2^q),
     * but can be calculated using addition and
     * bit shift operations.
     */
    unsigned change = (1<<(q+q))+ (s<<(q+1));
    if(a >= change) {
      a -= change;
      s |= 1<<q;
    }
    q--;
  }
  return s;
}

/***************************************
 * Verify my square root implementation
 * for all 32-bit unsigned integers
 **************************************/
int main(int argc, char *argv[])
{
  unsigned a=0;

  do {
    unsigned s = my_sqrt(a);
    if(s*s>a || s*s+2*s+1 < a) {
      printf("ERROR %u: %u\r\n",a, my_sqrt(a));
    }
    a++;
  } while(a != 0);
  return 0;
}


[slateraptor]

Warren actually notes this algorithm in his book. If you don't already own a copy, I think you'll def appreciate it from a practical utility perspective.

[Psi]

Here's an equally awesome inverse square root function which is super fast and will leave you scratching your head as to how it works,

Code: [Select]

float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;

x2 = number * 0.5F;
y  = number;
i  = * ( long * ) &y;                       // evil floating point bit level hacking
i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
y  = * ( float * ) &i;
y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
// y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed for faster but less accurate result

return y;
}

[coppice]

The SAR approach is quick if you have a fast multiplier, but in an FPGA implementation there are certainly denser ways to do it.

[legacy]

in the loop condition, it seems there is a missing char, probably "1"
while (a != 0); ----> while (a != 10);

Code: [Select]

int main(int argc, char *argv[])
{
    unsigned a = 0;

    a = 0;

    do
    {
        unsigned s = my_sqrt(a);
        printf("sqrt(%d)=%d\n", a, s);
        if (s * s > a || s * s + 2 * s + 1 < a)
        {
            printf("ERROR %u: %u\r\n", a, my_sqrt(a));
        }
        a++;
    } while (a != 10);
    return 0;
}


interesting method

[macboy]

Here's an equally awesome inverse square root function which is super fast and will leave you scratching your head as to how it works,

Code: [Select]

float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;

x2 = number * 0.5F;
y  = number;
i  = * ( long * ) &y;                       // evil floating point bit level hacking
i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
y  = * ( float * ) &i;
y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
// y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed for faster but less accurate result

return y;
}

The above comes from Quake III source, but that is not the origin of the algorithm. AFAIK, nobody has figured out exactly who came up with the algorithm and in particular the constant: 0x5f3759df.
One should note that the above code gives an estimation of the square root and would probably fail an IEEE754 validation. I would be curious to see how many NR iterations would be necessary to get the estimation close enough.

[f1rmb]

in the loop condition, it seems there is a missing char, probably "1"
while (a != 0); ----> while (a != 10);

Code: [Select]

int main(int argc, char *argv[])
{
    unsigned a = 0;

    a = 0;

    do
    {
        unsigned s = my_sqrt(a);
        printf("sqrt(%d)=%d\n", a, s);
        if (s * s > a || s * s + 2 * s + 1 < a)
        {
            printf("ERROR %u: %u\r\n", a, my_sqrt(a));
        }
        a++;
    } while (a != 10);
    return 0;
}


interesting method

Depending of the platform's size of unsigned, it's may be use to iterate for each bit, e.g on 8 bits, 255+1 == 0

Cheers.

Edit: just noticed the code,  it's from main(), so it's a burning test (iterates INT_MAX +1)

[miguelvp]

i  = 0x5f3759df - ( i >> 1 );   

Hex: 0x5f3759df is float: 13211836172961054720.0

Code: [Select]

Binary: 0101 1111 0011 0111 0101 1001 1101 1111

IEEE representation
sign exponent     mantissa
0    10111110 (1.)01101110101100111011111

Positive
2^(exponent-127)*[1.]mantissa

2^(190-127)*(1+1/4+1/8+1/32+1/64+1/128+1/512+1/2048+1/4096+1/32768+1/65536+1/131072+1/524288+1/1048576+1/2097152+1/4194304+1/8388608)

2^63*1.43243014812469482421875

9223372036854775808 * 1.43243014812469482421875

i >>1 is more complicated:
you are shifting the exponent to the right  2^128 ends up to 2^0 and adds 0.5 to the mantissa which is halved.

For example if the exponent was:
exp
128 -> 0,  1+0.5 + (mant/2)
127 -> 0, 1+mant/2
126 -> -1, 2^-1+0.5 + (mant/2)
125 -> -1, 2^-1 + (mant/2)
etc.

So it's doing something like:
i  = 13211836172961054720 - ((i*2^-128)/2+0.5*(i_exp%2))

but you have to extract the exponent from number, of course a negative number will have some trouble because the sign bit will shift into the exponent.

What is the significance of that formula? no idea, but it has to come from somewhere.

Edit: a quick google search for 13211836172961054720 gives this link where it explains it and has references including a paper that explains it further, it was attributed to Carmack but it goes back to 1974 PDP-11 UNIX.

http://mrob.com/pub/math/numbers-16.html

[ralphd]

Nice.  Anyone got an equally simple ln function or e^x?

[mikerj]

Nice.  Anyone got an equally simple ln function or e^x?

Log base 2 is pretty easy, the integer result is simply the bit number of the highest set bit in the word, and the rest of the word can be used to get an interpolated fractional result from a small table.  Once you have Log base two, just apply the standard log identity to change it to any other base.

2^x is the same, but in reverse.

[hamster_nz]

in the loop condition, it seems there is a missing char, probably "1"
while (a != 0); ----> while (a != 10);

I'm relying on 'a' wrapping from 0xFFFFFFFF to 0x0, allowing me to test all values.

Oh, and in another stupid-but-this-time-unintentional 'trick' I am relying on the optimizer to make the  "s * s + 2 * s + 1 < a" work, when really it shouldn't.

When s = 0xFFFF, then (s+1)*(s+1) = 0 and it was spitting errors, so I rewrote it as  s*s+2*s+1, which  should also equal 0, and it should also spit out errors during the test.

But the optimizer converts it into  "s * s + 2 * s  <= a", saving an addition, and as 0xFFFF * 0xFFFF + 2*0xFFFF = 0xFFFFFFFF it works by luck.

[chickenHeadKnob]

Nice.  Anyone got an equally simple ln function or e^x?

CORDIC versions of the elementary functions exist, ln and e^x can be indirectly calculated after computing inv TANH and SINH,COSH see Jack E Volder's original paper or these pdf's:

[url=http://casdc.ee.ncku.edu.tw/class/CA/CH22.pdf]http://casdc.ee.ncku.edu.tw/class/CA/CH22.pdf[/url]

archives.math.utk.edu/ICTCM/VOL11/C027/paper.pdf

edit: i futched up links

[Rasz]

   i  = 0x5f3759df - ( i >> 1 );               // what the fuck?

mm quake )
ze history on that one:
http://www.slideshare.net/maksym_zavershynskyi/fast-inverse-square-root

[coppice]

   i  = 0x5f3759df - ( i >> 1 );               // what the fuck?

mm quake )
ze history on that one:
http://www.slideshare.net/maksym_zavershynskyi/fast-inverse-square-root

If you want the history of the algorithm you need to go back much farther. I really don't know just how far, but Gauss is usually a good guess for who originally came up with things like that. :-)

[miguelvp]

If you just need the relation with another number and not the actual value, it's faster to square the other value and avoid doing the square root to begin with. Using distance squared is a common practice in computer graphics.

[nitro2k01]

If you want the history of the algorithm you need to go back much farther. I really don't know just how far, but Gauss is usually a good guess for who originally came up with things like that. :-)

Nope! Look closer. The real trick here is not the Newton-Raphson successive approximation which, as you mention, is centuries old. The real trick is the dirty casting from float to int to float, and doing an int operation on the float data. This requires a deeper understanding on how floats are stored bitwise, and saves one iteration in the successive approximation, which cuts the number of float multiplications in half.

[mikerj]

If you want the history of the algorithm you need to go back much farther. I really don't know just how far, but Gauss is usually a good guess for who originally came up with things like that. :-)

Something tells me the IEEE floating point format had yet to be invented 200 odd years ago

[ralphrmartin]

Here's an equally awesome inverse square root function which is super fast and will leave you scratching your head as to how it works,

Here's my awesome inverse square root function

function invsqroot(x)
return x*x
end

Did I miss something?

[Rasz]

Here's my awesome inverse square root function

function invsqroot(x)
return x*x
end

Did I miss something?

hi school?

[coppice]

If you want the history of the algorithm you need to go back much farther. I really don't know just how far, but Gauss is usually a good guess for who originally came up with things like that. :-)

Something tells me the IEEE floating point format had yet to be invented 200 odd years ago

IEEE754 only influences the bit pattern used. Not the principle.

[mikerj]

If you want the history of the algorithm you need to go back much farther. I really don't know just how far, but Gauss is usually a good guess for who originally came up with things like that. :-)

Something tells me the IEEE floating point format had yet to be invented 200 odd years ago

IEEE754 only influences the bit pattern used. Not the principle.

As someone else already mentioned, the novel part of this function is optimisation made by performing specific integer operations on a floating point value, rather than the actual Newton method of finding roots.

Discounting the IEEE format, where there any other binary floating point formats 200 years ago that would have enabled this to be done?

[nitro2k01]

IEEE754 only influences the bit pattern used. Not the principle.

No. Again, the optimization has to do with doing an integer operation on a float binary value. This single addition operation approximates, for the range of values used, the first approximation term well enough to be useful. This works because the float to int binary cast approximates (iirc) a log2 operation on the number. This is a non-trivial and novel optimization. Saying that this was known 200 years ago is like saying the Newton-Raphson method was known 3000 years ago just because the Babylonians could do addition and multiplication. After all, Newton-Raphson is "just" addition and multiplication, at least if you ignore the framework of calculus necessary to generalize the theory and come up with a representation for a given function.

[coppice]

If you want the history of the algorithm you need to go back much farther. I really don't know just how far, but Gauss is usually a good guess for who originally came up with things like that. :-)

Something tells me the IEEE floating point format had yet to be invented 200 odd years ago

IEEE754 only influences the bit pattern used. Not the principle.

As someone else already mentioned, the novel part of this function is optimisation made by performing specific integer operations on a floating point value, rather than the actual Newton method of finding roots.

Discounting the IEEE format, where there any other binary floating point formats 200 years ago that would have enabled this to be done?

The principle is to use log approximations. John Napier and others first used logs in the 16th century, and their value as shortcuts in all sorts of mathematical operations was quickly realised. A lot of log related tricks have long histories. That makes me strongly suspect someone will have done something like the software trick as a paper shortcut long ago, including their use to choose sensible starting approximations for interative algorithms. From the first floating point computers in the 1960s people were documenting all sorts of tricks to prime approximations based on the log nature that can be exposed by the interplay of integer and floating formats.

[nitro2k01]

That makes me strongly suspect someone will have done something like the software trick as a paper shortcut long ago, including their use to choose sensible starting approximations for interative algorithms. From the first floating point computers in the 1960s people were documenting all sorts of tricks to prime approximations based on the log nature that can be exposed by the interplay of integer and floating formats.

I'll still have to disagree somewhat. The main use for computers at that time was scientific/engineering/economics where a correct result generally trumped everything else. In a graphics renderer, you may not care about a 1%, 2% or even 5% error in extreme cases, because no one will complain. Especially if it means that the program can run at 60 FPS on someone's computer instead of 30 FPS. The idea of throwing accuracy out the window for raw speed, I argue, is younger than the '60s. In scientific calculations, you can never guarantee that the number won't be used in further calculations. In multimedia programming, you know that the display or speaker output is the final destination of the data, and the result is also not entirely critical. (Ie, the integrity of lives, buildings or economies don't depend on it.)

As for paper shortcuts, I doubt anyone would have used floating point tricks on paper, since binary floating is not a very natural notation system for humans. Of course, there are various other paper shortcuts.

[miguelvp]

The voodoo has nothing to do with  anything known to anyone in over 100 years ago.

The right shift not only does divide the mantissa in half but it adds 0.5 to it if the biased exponent is even it then multiplies the number by 2^-64

There was no mathematical representation that would allow doing all of that in the past.
Edit: plus there is an implied (1.) on the mantissa that is not represented on the shift, so it's (mantissa-1)/2+1

Then the subtraction goes further to subtract two numbers bit by bit with different biased exponents, so in theory what is doing is multiplying the number by 2^(numExp-ConstExp) doing the math, then shifting everything back and getting a good approximation that would cut on the iterations needed by using the newton method alone.

The constant is not arbitrary, I guess they started using log base 2 of 127, but then changed it to work with a specific range of numbers, seems the intention was to make it work with 55,300 in the middle but whoever did the program transposed the last two nibbles. But there is a number for which this approximation is tuned for, close enough to 55,300. I'll double check but I think it was 55,299.23 or something like that so it seems the intention wasn't targeted for that number. Edit: Correction, it waas 55299.32 I guess I have the same problem as Carmack

Yeah, if you plug 55,299.32 (0x47580352) that due to mantissa bits limitations is really 55299.32031, gives you the almost exact 1/sqrt(55299.32) = 0.00425246|12508714199066162109375, which is the closest representation in 32bit floats of the real answer (0.00425245870102170015628783736641)

The more you digress from that value (55299.32) the error increases, but by very little.

I'm surprised they didn't do the extra steps to use (n-55299.32) to get an extra compensation from the drift making it more accurate, and probably not needing the Newton method to bring the number more into convergence. Not a linear divergence (like 0x40 in the 65535 range and 0x20 in the 55300 range per integer increase) but seems like it would help get the solution closer with an extra multiplication so still no divisions.

Edit: What I meant by transposing the last two nibbles I mean the constant I think it was intended to be 5F3759FD which it will give the exact 1 over square root of 55300, maybe Carmack is a bit dyslectic.

[T3sl4co1l]

Nice.  Anyone got an equally simple ln function or e^x?

Yes.  Without using a table (potentially a large one, unwieldy for 8-bit MCUs), you can calculate the remaining bits by left-shifting the operand (N bits) until you have a 1.(N-1) fixed point number (the number of shifts is the integral part of the result N - Lg(operand)), then iterating the subtraction and squaring of the value.  Namely: each time you square the result, if the squaring generates a carry, then the result is greater than 2.  So we subtract 2 from the result (i.e., shift the carry into our output register as the next fractional bit of the log result) and continue.  After N passes, we could continue, but presumably, no more information can be had from it, so that the output has the same number of bits as the input.

Tim

[coppice]

That makes me strongly suspect someone will have done something like the software trick as a paper shortcut long ago, including their use to choose sensible starting approximations for interative algorithms. From the first floating point computers in the 1960s people were documenting all sorts of tricks to prime approximations based on the log nature that can be exposed by the interplay of integer and floating formats.

I'll still have to disagree somewhat. The main use for computers at that time was scientific/engineering/economics where a correct result generally trumped everything else. In a graphics renderer, you may not care about a 1%, 2% or even 5% error in extreme cases, because no one will complain. Especially if it means that the program can run at 60 FPS on someone's computer instead of 30 FPS. The idea of throwing accuracy out the window for raw speed, I argue, is younger than the '60s. In scientific calculations, you can never guarantee that the number won't be used in further calculations. In multimedia programming, you know that the display or speaker output is the final destination of the data, and the result is also not entirely critical. (Ie, the integrity of lives, buildings or economies don't depend on it.)

As for paper shortcuts, I doubt anyone would have used floating point tricks on paper, since binary floating is not a very natural notation system for humans. Of course, there are various other paper shortcuts.

There is no accuracy limitation. The graphics code is limited in accuracy, because they only run the iteration once or twice. The trick is to get a really good priming value for the loop, so it converges really quickly. That benefit is there for the crude graphics case, or with a few more iterations for the precision case.

You might want to read a little more about the kind of calculations computers did a century ago, when computer was a job title, and not a machine.

[ralphrmartin]

Here's my awesome inverse square root function

function invsqroot(x)
return x*x
end

Did I miss something?

hi school?

Inverse  square root is square. Maybe you missed something...

[miguelvp]

the negative value of the root degree makes that happen, not the inverse.

Edit: the term inverse in math is not always 1/ something in any event

[ralphrmartin]

the negative value of the root degree makes that happen, not the inverse.

Edit: the term inverse in math is not always 1/ something in any event

Inverse in advanced mathematics does not (always) mean 1/ something.

The inverse of a function f is a function, written f^-1, that undoes the effect of function f.

Trust me. I am a Chartered Mathematician as well as a Chartered Engineer. 

[miguelvp]

True 1/something as the inverse only applies to multiplication but all this is besides the point.

The function in question computes the (multiplication)inverse of the square root of a given number, so it's not x^2

Bringing semantics won't change that.

Edit: But when I first read your answer I thought it was obvious you said it in tongue in cheek.

[mikerj]

You might want to read a little more about the kind of calculations computers did a century ago, when computer was a job title, and not a machine.

You are still not seeing/understanding the principal operation that makes this function novel and fast, which was well described by Miguel above.  This is fundamentally a bit twiddling hack to reduce the number of iterations required, and is certainly not something you would do if you were manually using Newton Raphson or other methods for approximating roots.

[miguelvp]

I tracked down the original code (but using a different magic number) from:
Jim Blinn, Floating-point tricks, IEEE Comp. Graphics and Applications 17, no 4, 1997

Someone had a scan of the article but it's no longer available so Wayback Machine to the rescue:
https://web.archive.org/web/20130309073539/http://rufus.hackish.org/~rufus/FPtricks.pdf

Or get it from IEEE Xplore, but I haven't been an IEEE member for at least the last 15 years (I do have at least one paper published in there from 1994, however not in computer graphics, And I'm not sure I want to link it because it's really not relevant.) Anyways I'm not willing to pay $31 or even the $13 for members so I got it from the previous link.
http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=595279&queryText=jim+blinn&newsearch=true&searchField=Search_All

So we have to thanks Microsoft's Research for that Gem, even if there was an earlier but more obscure use
http://minnie.tuhs.org/cgi-bin/utree.pl?file=V5/usr/source/s3/sqrt.s (Dennis Ritchie 1974 as a support C library for the PDP-11)

Binn does attribute some ties to Graphic Gems 4 (which I do have but I donated all of them to work but I have access to them), but the actual Quake code as it's coded seems to have derived from that IEEE publication, and the magic number tweaked by Gary Tarolli (3DFX) probably to make it work better for the range they needed.

I will check on Monday to see if there is any similar code on Graphic Gems 4, I did look at Graphic Gems 1 on the past two days but didn't find anything relevant, but since Blinn cites volume 4 I bet there is something in there.

The constant on the IEEE paper is:
(OneAsInteger + (OneAsInteger>>1))

OneAsInteger is 0x3F800000 making the constant 0x5F400000 but as mentioned earlier that constant probably wasn't good enough for the range so it was tweaked.

For the meaning of: (OneAsInteger + (OneAsInteger>>1)) you have to read the paper and see what:
int i = (AsInteger(x)>>1)+(OneAsInteger>>1);

That gives you an approximation of the square of x, so it seems that adding 0x3F800000 as an interger performs the 1/x operation.

I guess I'll read that article in more detail because he (Jim Blinn) goes on to explain a lot of what is going on.

Edit: The article starts with this:

Well, I recently learned about some interesting properties
of the IEEE floating-point representation that sort
of mixes the two. To see (and believe) how this works,
we must first review this representation.

So the knowledge comes from somewhere else? Monday I'll check all of the Graphics Gems 1-4 (I have more but it has to be on the ones pre 1997)

[ralphrmartin]

OK, I'll stop joking around, and just note that the proper name for what someone called inverse square root is actually reciprocal square root.

[mikerj]

I tracked down the original code (but using a different magic number) from:
Jim Blinn, Floating-point tricks, IEEE Comp. Graphics and Applications 17, no 4, 1997

Someone had a scan of the article but it's no longer available so Wayback Machine to the rescue:
https://web.archive.org/web/20130309073539/http://rufus.hackish.org/~rufus/FPtricks.pdf

Thanks for posting that, it's a really interesting article.

[dom0]

If you want the history of the algorithm you need to go back much farther. I really don't know just how far, but Gauss is usually a good guess for who originally came up with things like that. :-)

Nope! Look closer. The real trick here is not the Newton-Raphson successive approximation which, as you mention, is centuries old. The real trick is the dirty casting from float to int to float, and doing an int operation on the float data. This requires a deeper understanding on how floats are stored bitwise, and saves one iteration in the successive approximation, which cuts the number of float multiplications in half.

I think the first publicly available occurrence of this specific method is somewhere in the original UNIX (for PDP 11) sources. Mid 70s.

[miguelvp]

I think the first publicly available occurrence of this specific method is somewhere in the original UNIX (for PDP 11) sources. Mid 70s.

As stated earlier but that used an unknown magic number stored at 0x4E84 and probably was done by Dennis Ritchie

...
So we have to thanks Microsoft's Research for that Gem, even if there was an earlier but more obscure use
http://minnie.tuhs.org/cgi-bin/utree.pl?file=V5/usr/source/s3/sqrt.s (Dennis Ritchie 1974 as a support C library for the PDP-11)
...