Truncated match.
PICList
Thread
'Floating point square root: Was: integer square ro'
1996\04\25@141855
by
John Halleck

My appologies to those that don't remember back far enough to remember
the thread that I promised to post this in.
This is the code to quickly compute sqrt(x*x + y*y)
It produces one significant digit in one interation,
5 significant digits in two iterations
20 significant digits in three interations,
62 significant digits in four iterations, etc.
The article (Referenced in the code) give versions that converge
much faster (for the same number of divides) for uses where division
is expensive (such as multiple precision arithmatic packages, and
PIC processors). But I don't have implementations for those around,
and I did have this around.
This is Fortran code from way back, but I assume that most folk here
can follow it. (Sorry, I've never implemented it on a PIC.)
At the time this came out, on all the 36 bit machines I had access to,
this routine was significantly faster than computing sqrt(x*x+y*y) in
the usual manner.
It's selling points were speed, and the lack of overflow unless the
answer overflowed. (I.E. if the answer fit on the machine, you got
it even if x*x+y*y might be larger than the machine could handle.)

FUNCTION ENORM2 (X, Y)
C
C This algorithm is from IBM J. Res. Development, Vol27 No 6, Nov 1983
C This code written by John Halleck April 1984
C
C Compute SQRT ( X*X + Y*Y )
C (Mathmaticaly this is the Euclidean vector norm for 2 dimensions)
C
C No computation in this routine will cause overflow unless the answer
C itself is too large to fit in the machine!
C
C *** Warning *** If the larger of X and Y is very close to
C floating point underflow the routine should return it instead of computing
C the answer. (This is a VERY rare case)
C
REAL X, Y
C The vector.
REAL GUESS, ERROR, DX, DY, R, S
C Current guess, residue, delta X, delta Y, and temporaries.
C
C Set up first guess.
DX = ABS (X)
DY = ABS (Y)
GUESS = AMAX1 (DX, DY)
ERROR = AMIN1 (DX, DY)
C
C The following code is the expanded form of a loop around
C R = (ERROR/GUESS)**2
C S = R / (4.0+R)
C GUESS = GUESS + 2.0 * S * P
C ERROR = ERROR * Q
C The termination condiction, if you don't unroll the loop, should
C be that R is numbericly insignificant with respect to 4.0
C
C 1 iteration gives 1 significant figure in the result.
C 2 iterations gives 5 significant figures,
C 3 iterations give 20 significant figures.
C 4 iterations gives 62 significant figures.
C
C The loop invariant is that SQRT(GUESS**2+ERROR**2) is
C SQRT (X**2+Y**2)
C
R = (ERROR/GUESS) ** 2
S = R / (4.0 + R)
GUESS = GUESS + 2.0 * S * GUESS
ERROR = S * ERROR
C
R = (ERROR/GUESS) ** 2
S = R / (4.0 + R)
GUESS = GUESS + 2.0 * S * GUESS
ERROR = S * ERROR
C
R = (ERROR/GUESS) ** 2
S = R / (4.0 + R)
GUESS = GUESS + 2.0 * S * GUESS
ERROR = S * ERROR
C
ENORM2 = GUESS
RETURN
END
More... (looser matching)
 Last day of these posts
 In 1996
, 1997 only
 Today
 New search...