Martin Guy, June 1985
Mr Woo's algorithm on a decimal abacus is a follows:
Put the number you want to root in the right hand end of the abacus. Count columns in twos from the decimal point, until you get to the greatest significant pair of columns. This will be the pair of columns you will start operating on first, as if the decimal point were just after the right hand column of the pair.Once the algorithm is reduced to a binary equivalent, dealing with pairs of bits at a time, nice properties drop out like:
Put 1 bead up on the leftmost column. This is worth 1 unit.
Subtract 1 unit from the MSP (Most Significant Pair of columns) as if they represented a number from 1 to 99.
Add two to the one on the left to make three.
If there are still three worth left in the MSP, subract three from it and add another two to the workspace to make five.
Carry on the cycle of subtracting the value of the workspace from the MSP and adding two to the workspace each time you could.
The last time you subtract the workspace from the MSP, just add one, not two.
Now multiply the workspace by ten and add one.
Move the MSP market two columns to the right.
When you have got down to the last two columns, if the original number was an exact square, the last subtraction should exactly wipe out the residue of the original number. If not, there is your remainder.
To get the root value, add a final one to the workspace, and divide it by two. There's your square root.
The loop subtracting the workspace from the original number can only ever do the loop contents once, or not at all, so it becomes a binary test.Tests of the integer version in compiled C on the VAX against the VAX floating point one in the new hand-crafted assembly language -lnm show the integer one to be six times as fast as the one in the library.
The bizarre adding and subtracting of one and two at various points can all be greatly reduced by using a floating decimal point in both registers.
Remember that incrementing a binary floating point exponent and a 1-bit shift can instantly align a floating point number to a 2-bit boundary. Change in the mantissa is simply achieved by multiplying it by two!
Here is the code in C: sqrt.c
Here, in Occam1: sqrt.occ
in Imp: sqrt.imp
in KRC: sqrt.krc
and in Miranda (courtesy of DAT): sqrt.m