算盘算法开根号

Square Root by Abacus Algorithm

Martin Guy, June 1985
after
C. Woo

ABSTRACT

A fast algorithm for calculating integer (and hence floating point) square roots is presented, with implementations in some popular and unpopular programming languages.

Mr Woo

It all started when I was given a fifteen-row abacus by my grandmother, with a book by a Mr C. Woo, a name almost too good to be true, on how to work it to do addition, subtraction, multiplication, division and square and cube roots. I was amazed. With some practice, I could do 10-digit square roots, giving remainder if not exact, in five or ten minutes.

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.

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.

Iterate.

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.

Once the algorithm is reduced to a binary equivalent, dealing with pairs of bits at a time, nice properties drop out like:

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.

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.

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.

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!

C代码：

/*
 *	Square root by abacus algorithm, Martin Guy @ UKC, June 1985.
 *	From a book on programming abaci by Mr C. Woo.
 *	Argument is a positive integer, as is result.
 *
 *	I have formally proved that on exit:
 *		   2		   2		   2
 *		res  <= x < (res+1)	and	res  + op == x
 *
 *	This is also nine times faster than the library routine (-lm).
 */

int
sqrt(x)
int x;
{
	/*
	 *	Logically, these are unsigned. We need the sign bit to test
	 *	whether (op - res - one) underflowed.
	 */
	
	register int op, res, one;

	op = x;
	res = 0;

	/* "one" starts at the highest power of four <= than the argument. */

	one = 1 << 30;	/* second-to-top bit set */
	while (one > op) one >>= 2;

	while (one != 0) {
		if (op >= res + one) {
			op = op - (res + one);
			res = res +  2 * one;
		}
		res /= 2;
		one /= 4;
	}
	return(res);
}

原始出处：http://medialab.freaknet.org/martin/src/sqrt/

更详细的在：

https://en.wikipedia.org/wiki/Methods_of_computing_square_roots

其中算盘算法x*(20p+x)是因为：

(10*p+x)*(10*p+x) - 10p*10p=100p*p+20px+x*x-100p*p=20px+x*x=x*(20p+x)