I’m having a hard time looking this up. Do any of you know a fast square root for 16:8 fixed point math? In C++ preferably, thanks.
If nobody knows of any, no problem.
Type in “fixed point square root algorithm” into Google (the phrase, not just the words), and click a little in the second result to find this.
Thank you. I tried several variations of “Fixed Point”, “Square Root”, and Algorithm, but never came across that link, or any like it. Thanks again.
Or wait, on second thought, just looking at it…I’ll give it a try…
Just what I thought, jibberish. Sorry.
Oops, turns out I need to take that second post back, it works fine after all.
It’s late, and I need to go to bed now, as evidence by my brain turning to goo.
Sorry about calling the link jibberish.
Actually, I get various degrees of success depending on the number. Oh well, maybe it’s me. Maybe I’ll just have to Google around a little bit more tomorrow.
You do realize that that’s an iterative algorithm, right? Each time you run it, you take an estimate for the square root, and get a better estimate out. If you then use that better estimate in the algorithm, you’ll get a better one still, and so forth. Eventually, you’ll get close enough to your answer. Really, the algorithm should also include some criterion for when to stop, but that’s easy enough to do in this case: When x and a/x are very close together (within your tolerances of error), that’s when you stop.
Probably not the most efficient way unless you will need arbitrary roots, but it does offer a fixed processing time:
Taking the base-2 log. of fixed point binary is pretty simple.
Gak, hit submit too soon.
OK,
Base 2 logs work nice for binarys in the same way that base 10 logs work well for
deca-digit handicapped humans.
So the base two log can be found by just counting down from MSB untill you find the first one bit, which gets you the exponent…OK, so that part isn’t fixed timing, but it is a trivial loop exicuted at most 15 times.
Then you use a look up table to find the mantissa (terminology?). The table doesn’t need a lot of entries…16 is probably enough, and you may need to then interpolate. If the exponent was less than 5, then you don’t need to interpolate.
Divide the log by two (bit shift), and take the anti-log, and you’ve got your root. If you only need square roots, you can pre shift the table entries and save a step. If you divide by n, then you get the nth root instead of a square root.
IIRC, Intel used to have an app note on this in thier 8-bit microcontrollers databook…back when giants walked the earth.
I had one that I came up with which you can try against the one linked to.
Simply, you create an array where the item at each position is the square of the index. Then when you go for the root, you just do a binary search for the value in the array and return the index it was found at.
Note that, to save memory, I only allowed the root of 20-bit values, so the array length was only 2[sup]10[/sup] big.
Also note that not every value is in the array, so you need to do a search for the **nearest/b] value.
Anyways, my usenet thread can be found here with the source.
Last night I was extremely tired, and my brain turned to mush, as should be evident by my incoherent ramblings. I wonder if that’s what being drunk feels like. 
Anyway, I thought I had a handle on the algorithm, but when I tested it, I apparently had done something wrong. Now that I’m actually awake, I can study it some more to see what I need to correct.
Thank you, I’ll check it out.
This algorithm is actually equivalent to Newton’s method applied to the function f(x) = x[sup]2[/sup] - a, so yes, it works, and pretty fast. On the other hand, not every initial guess will converge to the solution. Maybe that’s why it doesn’t work all the time for you? Try taking a value close enough to the square root of a for an initial guess.
Yeah. Actually, thinking about it, I don’t think that kind of algorithm isn’t what I’m looking for. Well, to be honest, I don’t actually need it right now, but in the future I’m sure I’ll have use for one. The problem whith that algorithm is, the guessing. The way I’ll be using a square root, it’ll all be variables, which means I won’t know what exact value I’ll be using, which means I won’t know that to guess.
Well, every iterative algorithm to solve your problem will require you to start with an initial guess. In fact, with this kind of quadratic polynomial, I think that every initial guess greater than 0 will eventually converge to the solution. So it’s not too bad.
For more information on Newton’s (Newton-Raphson’s) method:
From Wikipedia
From MathWorld
Ah ha, I found my fatal flaw.
At first I was having the algorithm repeat itself while x * x != a. It worked for some numbers, but no others.
Then I tried having it repeat itself while x * x < a. Didn’t work either.
But having it repete while x * x > a, that works perfectly.
As Chronos said, what is usually done is to fix a given error tolerance (you choose which one) and break the loop at iteration n when |x[sub]n[/sub][sup]2[/sup] - a| (or |x[sub]n[/sub] - a/x[sub]n[/sub]|, as he suggested) is smaller than this tolerance. I’m not sure why your stopping criterion would work. I’ll do an example: let’s say you want to find the square root of 55, and start with an initial guess of x[sub]0[/sub] = 7.4 (since the square root of 49 is 7 and the square root of 64 is 8 and 55 is somewhere between these numbers). Now, 7.4[sup]2[/sup] = 54.76, so with your stopping criterion, the algorithm would stop right away, but you’re still pretty far. Instead, take, for example, 10[sup]-6[/sup] as a tolerance. This is what you get (define f(x) = x[sup]2[/sup] - 55):
x[sub]0[/sub] = 7.4; |f(x[sub]0[/sub])| = 0.24
x[sub]1[/sub] = (x[sub]0[/sub] + 55/x[sub]0[/sub])/2 = (7.4 + 55/7.4)/2 = 7.41622; |f(x[sub]1[/sub])| = 2.63 x 10[sup]-4[/sup]
x[sub]2[/sub] = (x[sub]1[/sub] + 55/x[sub]1[/sub])/2 = (7.41622 + 55/7.41622)/2 = 7.41620; |f(x[sub]2[/sub])| = 3.14 x 10[sup]-10[/sup]
and at this point the algorithm stops. The square root of 55 is therefore close enough to 7.41620. (I did not include all the decimals I used, and there are ways to know up to what point your decimals are correct, such as doing an extra iteration and seeing which decimals don’t change anymore.)
Do you see how this works?
Note that in this particular case (by this I mean with this particular function f(x) = x[sup]2[/sup] - a), almost all your values of x[sub]n[/sub] will be larger than the exact solution. If you start with an initial guess x[sub]0[/sub] that is larger than sqrt(a), all the next iterations will also be larger than sqrt(a), while if you start with a x[sub]0[/sub] smaller than sqrt(a), then all the following iterations will be larger than sqrt(a). This can be easily proven and illustrated using the definition of Newton’s method and the graph of the parabola f(x) = x[sup]2[/sup] - a. I guess this is what you meant by using x[sup]2[/sup] > a as a “repeat” criterion. You start with an initial guess larger than sqrt(a) and repeat until the distance between x[sub]n[/sub][sup]2[/sup] and a is too small for your machine. But I’ve given an example of a case where it may fail, and it doesn’t allow you to choose your own tolerance for the stopping criterion.
Although my test was for floats, I’m going to actually use it on 16:8 fixed point. The smallest unit for such a fixed point number is .01 so, my initial guess is allways .01 because the way I’ll be using the algorithm, I’ll never know what the input will be, thus, I can’t guess on the unknown. I doubt I’ll be using actual numbers, but variables instead, which could be anything, how do I guess that? So always starting with the smallest number allowed (smallest positive number) is the safest way to go. And the initial calculation always gives a bigger result the number who’s square root you’re looking up, and each time the algorithm is repeated, the numbers get smaller and smaller, till they either match the number you’re looking up, or are slightly smaller.
Here’s my code for testing the algorithm. Works every time.
#include<stdio.h>
float A, B, SqrRoot(float C);
unsigned short D = 0;
main()
{
printf("%G", SqrRoot(.5));
printf(" Took - %u tries", D);
return 0;
}
float SqrRoot(float C)
{
float x = .01;
do{
x = (x + C/x)/2;
D++;
}while((x * x) > C);
return x;
}