Author Topic: Calculating square roots.  (Read 13170 times)

0 Members and 1 Guest are viewing this topic.

Offline hamster_nzTopic starter

  • Super Contributor
  • ***
  • Posts: 2803
  • Country: nz
Calculating square roots.
« on: June 29, 2015, 10:49:12 am »
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;
}
« Last Edit: June 29, 2015, 10:54:04 am by hamster_nz »
Gaze not into the abyss, lest you become recognized as an abyss domain expert, and they expect you keep gazing into the damn thing.
 

Offline slateraptor

  • Frequent Contributor
  • **
  • Posts: 833
  • Country: us
Re: Calculating square roots.
« Reply #1 on: June 29, 2015, 11:27:34 am »
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.
 

Online Psi

  • Super Contributor
  • ***
  • Posts: 9888
  • Country: nz
Re: Calculating square roots.
« Reply #2 on: June 29, 2015, 11:34:30 am »
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;
}
« Last Edit: June 29, 2015, 11:36:25 am by Psi »
Greek letter 'Psi' (not Pounds per Square Inch)
 

Offline coppice

  • Super Contributor
  • ***
  • Posts: 8605
  • Country: gb
Re: Calculating square roots.
« Reply #3 on: June 29, 2015, 12:31:01 pm »
The SAR approach is quick if you have a fast multiplier, but in an FPGA implementation there are certainly denser ways to do it.
 

Offline legacy

  • Super Contributor
  • ***
  • !
  • Posts: 4415
  • Country: ch
Re: Calculating square roots.
« Reply #4 on: June 29, 2015, 01:04:07 pm »
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
 

Offline macboy

  • Super Contributor
  • ***
  • Posts: 2252
  • Country: ca
Re: Calculating square roots.
« Reply #5 on: June 29, 2015, 01:39:33 pm »
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.
 

Offline f1rmb

  • Regular Contributor
  • *
  • Posts: 180
  • Country: fr
Re: Calculating square roots.
« Reply #6 on: June 29, 2015, 01:55:55 pm »
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)
« Last Edit: June 29, 2015, 04:17:21 pm by f1rmb »
 

Offline miguelvp

  • Super Contributor
  • ***
  • Posts: 5550
  • Country: us
Re: Calculating square roots.
« Reply #7 on: June 29, 2015, 05:20:20 pm »
Quote
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

« Last Edit: June 29, 2015, 05:47:51 pm by miguelvp »
 

Offline ralphd

  • Frequent Contributor
  • **
  • Posts: 445
  • Country: ca
    • Nerd Ralph
Re: Calculating square roots.
« Reply #8 on: June 29, 2015, 06:40:48 pm »
Nice.  Anyone got an equally simple ln function or e^x?
Unthinking respect for authority is the greatest enemy of truth. Einstein
 

Online mikerj

  • Super Contributor
  • ***
  • Posts: 3233
  • Country: gb
Re: Calculating square roots.
« Reply #9 on: June 29, 2015, 08:30:25 pm »
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.
 

Offline hamster_nzTopic starter

  • Super Contributor
  • ***
  • Posts: 2803
  • Country: nz
Re: Calculating square roots.
« Reply #10 on: June 29, 2015, 09:34:09 pm »
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. :D


Gaze not into the abyss, lest you become recognized as an abyss domain expert, and they expect you keep gazing into the damn thing.
 

Offline chickenHeadKnob

  • Super Contributor
  • ***
  • Posts: 1054
  • Country: ca
Re: Calculating square roots.
« Reply #11 on: June 30, 2015, 12:19:17 am »
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
« Last Edit: June 30, 2015, 12:42:08 am by chickenHeadKnob »
 

Offline Rasz

  • Super Contributor
  • ***
  • Posts: 2616
  • Country: 00
    • My random blog.
Re: Calculating square roots.
« Reply #12 on: June 30, 2015, 09:21:56 am »
   i  = 0x5f3759df - ( i >> 1 );               // what the fuck?

mm quake :))
ze history on that one:
http://www.slideshare.net/maksym_zavershynskyi/fast-inverse-square-root
Who logs in to gdm? Not I, said the duck.
My fireplace is on fire, but in all the wrong places.
 

Offline coppice

  • Super Contributor
  • ***
  • Posts: 8605
  • Country: gb
Re: Calculating square roots.
« Reply #13 on: June 30, 2015, 09:33:57 am »
   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. :-)
 

Offline miguelvp

  • Super Contributor
  • ***
  • Posts: 5550
  • Country: us
Re: Calculating square roots.
« Reply #14 on: June 30, 2015, 03:15:54 pm »
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.
 

Offline nitro2k01

  • Frequent Contributor
  • **
  • Posts: 843
  • Country: 00
Re: Calculating square roots.
« Reply #15 on: July 01, 2015, 11:29:23 am »
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. :)
Whoa! How the hell did Dave know that Bob is my uncle? Amazing!
 

Online mikerj

  • Super Contributor
  • ***
  • Posts: 3233
  • Country: gb
Re: Calculating square roots.
« Reply #16 on: July 01, 2015, 11:52:37 am »
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 :)
 

Offline ralphrmartin

  • Frequent Contributor
  • **
  • Posts: 479
  • Country: gb
    • Me
Re: Calculating square roots.
« Reply #17 on: July 02, 2015, 10:16:42 am »
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? :-DD
 

Offline Rasz

  • Super Contributor
  • ***
  • Posts: 2616
  • Country: 00
    • My random blog.
Re: Calculating square roots.
« Reply #18 on: July 02, 2015, 10:19:57 am »
Here's my awesome inverse square root function

function invsqroot(x)
return x*x
end

Did I miss something? :-DD

hi school?
Who logs in to gdm? Not I, said the duck.
My fireplace is on fire, but in all the wrong places.
 

Offline coppice

  • Super Contributor
  • ***
  • Posts: 8605
  • Country: gb
Re: Calculating square roots.
« Reply #19 on: July 02, 2015, 10:42:12 am »
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.
 

Online mikerj

  • Super Contributor
  • ***
  • Posts: 3233
  • Country: gb
Re: Calculating square roots.
« Reply #20 on: July 02, 2015, 11:40:34 am »
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?
 

Offline nitro2k01

  • Frequent Contributor
  • **
  • Posts: 843
  • Country: 00
Re: Calculating square roots.
« Reply #21 on: July 02, 2015, 12:24:39 pm »
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.
« Last Edit: July 02, 2015, 12:27:49 pm by nitro2k01 »
Whoa! How the hell did Dave know that Bob is my uncle? Amazing!
 

Offline coppice

  • Super Contributor
  • ***
  • Posts: 8605
  • Country: gb
Re: Calculating square roots.
« Reply #22 on: July 02, 2015, 12:58:35 pm »
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.
 

Offline nitro2k01

  • Frequent Contributor
  • **
  • Posts: 843
  • Country: 00
Re: Calculating square roots.
« Reply #23 on: July 02, 2015, 01:39:25 pm »
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.
Whoa! How the hell did Dave know that Bob is my uncle? Amazing!
 

Offline miguelvp

  • Super Contributor
  • ***
  • Posts: 5550
  • Country: us
Re: Calculating square roots.
« Reply #24 on: July 02, 2015, 02:27:01 pm »
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.

« Last Edit: July 03, 2015, 07:22:38 pm by miguelvp »
 


Share me

Digg  Facebook  SlashDot  Delicious  Technorati  Twitter  Google  Yahoo
Smf