The other day, I was playing around with fixed point maths - and for some reason I very much wanted to pull the square root out of a number. As I googled, I found pretty much no information. So I checked out the source code of Allegro, to see how they handled it in their fixed point math routines. They used two apporaches: On x86 they made a lookup for the first 8 bits, and shifted the result to the right position. This is fast, but gives you 2.25 significant numbers - so it's a pretty dirty solution. Besides, they used some x86-only assembly. So, how dothey do it on other platforms? They cast into float, and call sqrt(). This is where I called my father (a retired maths teacher), to find a better algorithm.
This is a method for solving square roots by hand (as they did it before they calculator). It's not harder than division, actually - it's a lot like division.
I'll start with a pretty big number: 123456789 Then we group the figures, with two in each group: 01 23 45 67 89 As you see I added a leading zero, this is because we need to get the decimal point between two groups.
Now we add two columns, one two the left of our number, and one to the right. The one to the left will add stuff up, and the one to the right will contain the result.
1 01 23 45 67 89 1 +1 - 1 2 0
We put the same figure in the left and right columns, and that must be the number which multiplied by itself is no bigger than the number in the first group (01 in our case). We also subtract that from the first group.
1 01 23 45 67 89 1 +1 - 1 21 0 23 1 + 1 - 21 22 2
In the next row, we shall move the second group down from our number, and hang it onto the difference in the middle column (0 in our case). Then we must find a number to put at the end of the sum in the left column. This number multiplied by the number we got when we put it at the end of the sum, shall be no bigger than the middle number (start with 1, and work yourself upwards).
1 01 23 45 67 89 1 +1 - 1 21 0 23 1 + 1 - 21 221 2 45 1 + 1 - 221 222 24
In the 5th row, there's another 1 - since 2*222 is way larger than 245, but 1*221 fits nicely. It seems I picked a boring number. ;)
1 01 23 45 67 89 1 +1 - 1 21 0 23 1 + 1 - 21 221 2 45 1 + 1 - 221 2221 24 67 1 + 1 - 2221 2222 246
And in the 7th row - we have yet another 1! I would never have imagined...
1 01 23 45 67 89 1 +1 - 1 21 0 23 1 + 1 - 21 221 2 45 1 + 1 - 221 2221 24 67 1 + 1 - 2221 22221 246 89 1 + 1 - 22221 22222 2468
And in the 9th row, we have 1 again. But now we are out of numbers - what now? Well, there is an infinite amount of 0 after the decimal point...
1 01 23 45 67 89 1 +1 - 1 21 0 23 1 + 1 - 21 221 2 45 1 + 1 - 221 2221 24 67 1 + 1 - 2221 22221 246 89 1 + 1 - 22221 . 222221 246800 1 + 1 -222221 2222221 2457900 1 + 1 -2222221 22222221 23567900 1 + 1 -22222221 222222220 134577900 0 + 0 - 0 2222222206 13457790000 6 + 6 -13333333236 2222222212 124456764
And so on, and so forth. If I do another line, my old casio won't be able to come along and check the numbers. But hopefully you get the idea how it's done manually by now. The result turns out to be 11111.11106 Do note that the number of the length of the integer part of the answer is as long as there are groups in the original number.
If we are to do this to fixed point numbers, we want to do it to binary numbers - not decimal. Is that a problem? No, as it turns out - it is even easier. When we do it with binary numbers, we never have to multiply the number from the left column before subtracting it from the middle column. Or at least, we always multiply by 1 or 0, which are no-ops.
In this function, I take a 64-bit argument and return a 32-bit result. The format of the 64-bit fixed is 32.32, and that gives us a 16.16 result.
typedef int32_t fixed_t;
typedef int64_t fixed64_t
fixed_t sqrt(const fixed64_t x)
{
// Can't pull a root out of a negative value
assert(x >= 0);
fixed_t result=0;
int group, sum=0, diff=0; // sum is the left column, and diff is the middle column
for (group=31; group>=0; --group) // The 64bit number is hadnled as 32 groups, group is the counter
{
// Get a new test sum, by shifting the old one left, and add one
sum=sum<<1;
sum+=1;
// Get a new diff,
// by shifting the old one left two steps,
diff=diff<<2;
// and add the current cell (masked with binary 11 == 3)
diff+=(x>>(group*2))&3;
// Check if the test sum is larger than the difference
if (sum>diff)
{
// Difference is not changed, and this bit is a 0
// in the result. We just correct the sum.
sum-=1;
}
else
{
// Set this bit,
result+=(1<<group);
// and calculate new vars
diff-=sum;
sum+=1;
}
}
return result;
}
You could possibly speed this function up by not starting at group 31, but at the first group with a bit set. But honestly - this is a twelve-liner. It does 32 iterations. Can it be done faster by introducing more complexity? OK, if you start with a 32-bit fixed point number, you only need to do 16 iterations. But not even that is much of a win.