-
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 usersCode: [Select]
/*****************************************
* my_sqrt.c: Implmenting an integer square
* root using only bit-shifts
*
* Author: Mike Field <[email protected]>
*
*****************************************/
#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;
} -
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.
-
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 SAR approach is quick if you have a fast multiplier, but in an FPGA implementation there are certainly denser ways to do it.
-
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 -
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,
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.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;
}
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. -
in the loop condition, it seems there is a missing char, probably "1"
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
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
Cheers.
Edit: just noticed the code, it's from main(), so it's a burning test (iterates INT_MAX +1) -
Quote
i = 0x5f3759df - ( i >> 1 );
Hex: 0x5f3759df is float: 13211836172961054720.0Code: [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
-
Nice. Anyone got an equally simple ln function or e^x?
-
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. -
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.
-
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 -
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. :-)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 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.
-
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.
-
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
-
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?
-
Here's my awesome inverse square root function
function invsqroot(x)
return x*x
end
Did I miss something?
hi school? -
IEEE754 only influences the bit pattern used. Not the principle.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.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
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? -
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. -
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.
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.
IEEE754 only influences the bit pattern used. Not the principle.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
Discounting the IEEE format, where there any other binary floating point formats 200 years ago that would have enabled this to be done? -
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. -
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.