So you need an integer SQRT function, but which to choose?

I've been analysing the performance of many different implementations of 16 bit integer square root, testing them exhaustively over all possible inputs, cycle counting them, then graphing the results. There are big differences in both performance and memory use:

https://github.com/TobyLobster/sqrt_test

I've omitted routines that calculate squares by adding successive odd numbers, since these are *way* too slow for anything but small numbers.

_________________
x86?  We ain't got no x86.  We don't NEED no stinking x86!

_________________
http://WilsonMinesCo.com/ lots of 6502 resources
The "second front page" is http://wilsonminesco.com/links.html .
What's an additional VIA among friends, anyhow?

Thanks for posting this work over here TobyLobster - I enjoyed the thread over on stardot!

Interesting results, for sure.

Michael J. Mahon wrote a very fast square root routine to replace the small but very slow exp(log(x)/2) Applesoft version. It performs its magic in 40-bit Microsoft floating point, but I have a crazy notion to simplify it way down for integers. I don't understand how it works, which is a bit of an impediment, and it's quite possible that its algorithm is already in one of the routines listed.

http://michaeljmahon.com/USR.SQR.html

_________________
Got a kilobyte lying fallow in your 65xx's memory map? Sprinkle some VTL02C on it and see how it grows on you!

Mike B. (about me) (learning how to github)

That sounds intriguing, I hadn't thought to convert a floating point SQRT to integer, but it sounds like it should work. It's not obvious to me if it will turn out to be similar to one I've already tested or not, but it will be interesting to find out! Let us know how you get on.

I am sure that many of the algorithms you've linked to in this thread are much more efficient than the algorithm I used to simultaneously calculate sqrt(S) and 1 / sqrt(S) for a control application that I implemented with a microprogrammed FPGA. It was based on a floating point algorithm, and I had to derive the integer form of the algorithm for my particular application where I specifically needed both the sqrt(S) and 1 / S. I derived the inverse of S by squaring the 1 / sqrt(S) term that the algorithm provided for free. I've described the derivation of my algorithm here.

_________________
Michael A.

Thanks for that article, Michael. I think I'll have to read that more than a couple of times to get the hang of the math.

In other news, I was playing around with MJM's MS float routine, and I found what I thought was an interesting little resting point:
Load the 16-bit integer into ACC, and it appears to come back as the square root, in 8.8 fixed point format! I still don't understand how it works, but it seems I can get just the truncated integer result by dumbing the ldx #16 down to ldx #8.

There appears to be room for improvement if the truncated result is all that's needed, but I'm chuffed at the current state of affairs. Can anyone with a better brain than mine try to explain what's going on here?

_________________
Got a kilobyte lying fallow in your 65xx's memory map? Sprinkle some VTL02C on it and see how it grows on you!

Mike B. (about me) (learning how to github)

Here's an attempt at an explanation of the algorithm.
Assume the input x is a fraction with 0 <= x < 1
In binary, this means the binary point is immediately to the left of the input value
Each iteration generates a single result bit
The first iteration determines the 2^-1 bit, the second the 2^-2 bit etc.
After n iterations, let y_n be the partial result
For iteration n+1 we need to decide the value of the 2^-(n+1) bit
The question which needs to be answered is as follows:

Is (y_n + 2^-(n+1))^2 >= x ?

If the answer is yes, the next bit is 1, otherwise it is 0

(y_n + 2^-(n+1))^2 = y_n^2 + 2^-n y_n + 2^-(2n+2)

This is a multiple of 2^-(2n+2), so at iteration (n+1), 2(n+1) fraction bits are involved
in the comparison, which means you need to take the next 2 input bits at each iteration

We rearrange the equations in terms of the error after n iterations:

e_n+1 = x - (y_n + 2^-(n+1))^2 = (x - y_n^2) - (2^-n y_n + 2^-(2n+2))
= e_n - (2^-n y_n + 2^-(2n+2))

Multiplying both sides by 2^2n gives

2^2n e_n+1 = 2^2n e_n - 2^n y_n - 2^-2

The terms on the RHS are as follows:

2^2n e_n = current error term shifted left so it is an integer, along with 2 extra input bits shifted into 2^-1 and 2^-2 positions
2^n y_n = current partial result shifted left so it is an integer
-2^-2 = an extra 1/4, which is where the CMP #$40 comes in

If the subtraction gives a borrow, the next result bit is 0 so we shift the partial result left 1
If the subtraction gives no borrow, the next result bit is 1 so we shift the partial result left 1 and add 1,
and update the current error value according to the result of the subtraction

Note that after n iterations, the difference between y_n and the true square root is guaranteed to be <2^-n
Let the true square root be r.

(r - y_n)(r + y_n) = r^2 - y_n^2
= x - y_n^2

e_n = (x - y_n^2) < 2^-n . 2 = 2^-(n-1)

So the error accumulator e_n grows by only 1 bit per iteration rather than 2

There's a similar algorithm for calculating cube roots. Here's a C implementation for 32 bit values:

@kernelthread: Thank you so much for taking the time and effort to write that up, and for solving the mystery of the #$40 constant! My 19-year-old brain would have gobbled up the theory, but my 56-year-old brain is struggling mightily.

Neat how the 6502's SBC instruction seems tailor-made for that, huh?

_________________
Got a kilobyte lying fallow in your 65xx's memory map? Sprinkle some VTL02C on it and see how it grows on you!

Mike B. (about me) (learning how to github)

I can't believe how close the SQRT routine is to WOZ's MOD routine.

This is from WOZ's Integer Basic. Note: In my assembler ":" is always a forward branch and "]" is always a backwards branch

Here is my implementation that is using 2 tables of 256 bytes and binary search. The code is in MADS assembler syntax for table generation.

Interesting approach I've not seen elsewhere. A couple of minor fixes are needed: You need "ldx r" before the rts in stop: to round the result down properly e.g. for SQRT(5), and inputs over 255*255 need handling properly (either with a special check at the start, such as:



or by adding BEQ stop into lower:



Timing shown in the graph below as sqrt11, in red

Thanks, very good points about the rounding. I needed it for an Euclidean distance and rounding problem was negligible for me, so I didn't catch it before.
The code is on the MIT license so feel free to use it or add it to your benchmark.