[net.math] need fast algorithm for reciprocal

geoff@burl.UUCP (geoff) (06/12/84)

I am working on a very primitive machine which lacks a division
operation.  What I really need is a reciprocal, so either type
algorithm would be fine.  I have shift, +,-, and *, as well as
logical operators, available.  I recall that IBM had an algorithm
many years ago which used multiplication to divide (sort of a binary
search -- multiply guess * number, compare to 1, adjust guess), but
I don't remember it well enough (my attempt at such an algorithm
required an avg of 20 iterations for reciprocals of numbers between 2
and 8000).

Any help from the net would be sincerely appreciated.

	thanks,
		geoff sherwood
		burl!geoff

	<"take care of the luxuries
		and the necessities take care
			of themselves">

ken@turtlevax.UUCP (Ken Turkowski) (06/20/84)

Suppose we want to calculate:			C = A / B.
This can be equivalently calculated as:		C = A * (1/B)
Consider the equation:				F(X) = (1/X) - B
The solution of:				F(X) = 0
is the reciprocal of B.

Use Newton's method to find the root:	X[i+1] = Xi - F(Xi)/F'(Xi)
					       = Xi - (1/Xi - B) / -Xi^2
					       = Xi(2 - B * Xi)

One can find the initial iterate, X0, from table lookup.  Convergence
is quadratic: i.e. the precision doubles with every iteration.  For
example, suppose the table lookup has 8 bits of accuracy; the first
iteration will produce 16 bits of accuracy, the second will produce 32.

The result will converge as long as the initial iterate, X0, is in the
range:

	0 < X0 < 2/B,		B > 0
	2/B < X0 < 0,		B < 0

-- 
Ken Turkowski @ CADLINC, Palo Alto, CA
UUCP: {amd70,decwrl,flairvax}!turtlevax!ken

johna@haddock.UUCP (06/26/84)

#R:burl:-48000:haddock:13000001:000:902
haddock!johna    Jun 25 16:16:00 1984

The following is a division routine to divide 2 floating point numbers
that have up to 31 bit mantissas.  It uses only shifts, ORs, ANDs and
subtraction.

#define LO31BITS 0x7fffffff
#define BIT0 0x1
divide(dividend, divisor, quotient,
	 sign_dvd, sign_dvi, sign_quo, exp_dvd, exp_dvi, exp_quo)

	long dividend, divisor, quotient, sign_dvd, sign_dvi, sign_quo;
	long exp_dvd, exp_dvi, exp_quo;
{
	register short i;

	sign_quo = sign_dvd ^ sign_dvi;
	exp_quo  = exp_dvd - expdvi;
	quotient = 0;
	for (i=0; i<31; i++)
	{       quotient = (quotient << 1) & LO31BITS;          /* shift up */
		if (dividend >= divisor)                /* can we subtract? */
		{       dividend -= divisor;                  /* then do it */
			quotient |= BIT0;            /* set lo-order bit on */
		}
		dividend = (dividend << 1) & LO31BITS;          /* shift up */
	}
}

Hope this helps.
		      ...!cca!ima!haddock!johna