| From: ROLL::USENET "USENET Newsgroup Distributor" 14-JUN-1984 22:14
To: HARE::STAN
Subj: USENET net.math newsgroup articles
Newsgroups: net.math
Path: decwrl!decvax!cca!g-rh
Subject: Calculation of reciprocal without division
Posted: Wed Jun 13 00:02:33 1984
{
A recent submission asked for a good algorithm for division,
given shift, multiply, and add instructions. The following is an
efficient and simple algorithm:
Suppose we have an approximation x to 1/N. Then we have
xN = 1+e which can be written as 1/N = x/(1+e). We have
1/(1+e) = 1-e+e^2-e^3+e-4 ...
= (1-e)(1+e^2)(1+e^4)(1+e^8)...
An estimate for x which is correct to the first bit can be found
by shifting N. This should be done so that xN lies in the range
[1/2,1]. The magnitude of the error term, e, is then bounded by
1/2. Four terms in the product are required for 16 bit accuracy
and five for 32 bit accuracy. The powers are formed by successive
squaring. A number of refinements are possible; the best algorithm
for a particular machine will depend on the exact instruction set
available and their execution times.
}
|
| From: ROLL::USENET "USENET Newsgroup Distributor" 21-JUN-1984 22:04
To: HARE::STAN
Subj: USENET net.math newsgroup articles
Newsgroups: net.math
Path: decwrl!turtlevax!ken
Subject: Re: need fast algorithm for reciprocal
Posted: Wed Jun 20 11:36:15 1984
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
|
| From: ROLL::USENET "USENET Newsgroup Distributor" 27-JUN-1984 22:09
To: HARE::STAN
Subj: USENET net.math newsgroup articles
Newsgroups: net.math
Path: decwrl!decvax!cca!ima!haddock!johna
Subject: Re: need fast algorithm for reciprocal - (nf)
Posted: Mon Jun 25 20:37:50 1984
#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
|
|
AN ALGORITHM FOR EXTENDED PRECISION DIVISION
Lynn Yarbrough, Development Methods
In working on projects on the PDP-11 which require arithmetic processing
in excess of that provided by the -11 hardware, I have had occasion to
develop some extended-precision capabilities. In particular, the
problem of computing the quotient and remainder when one
(Extended-Precision, or EP) integer is divided by another is a
challenge.
The following recursive algorithm for integer division has some pleasant
properties: based on primitive routines for EP addition, subtraction,
and comparison, the algorithm is quite short and reasonably fast (its
speed is proportional to the length, L, of the arguments and the speed
of the other subroutines). A surprising feature of the algorithm is
that it does not require the primitive operation of EP shifting.
The algorithm is mathematically based on the simple identity
A/B = 2*(A/2B)
in which the multiplications on the right-hand side can be computed by
addition. The parenthesized term will, after at most L repetitions of
the doublings, become less than one; at that point, the quotient is
zero and the remainder is the dividend. In unwinding the recursion, 1's
are injected into the quotient as required and the remainder is
correspondingly reduced. In the expression of the algorithm on the
following page, the operations :=, +, -, and the relational operators
are assumed to be applied to integers in EP representation.
Page 2
PROCEDURE LONGDIV (A, B, QUO, REM):
! (Inputs: A, B; Outputs: QUO, REM.)
BEGIN
LOCAL QQ, RR;
IF A < B
THEN ! (This always ends the recursion...)
BEGIN
QUO := 0; ! EP constant Zero
REM := A;
END
ELSE ! (The two EP doublings follow:)
BEGIN
LONGDIV (A, B+B, QQ, RR);
QUO := QQ+QQ;
REM := RR;
IF B <= RR
THEN ! The remainder is at most B too large and
! the quotient is one too small.
BEGIN
QUO := QUO + 1; ! EP constant One.
REM := REM - B;
END;
END;
END.
|