ComplexImpl.mesa
Last edited by Michael Plass, 29-Sep-81 9:14:51
Written by Michael Plass, 29-Sep-81
DIRECTORY Ieee USING [SingleReal], Vector USING [Mul], Complex;
ComplexImpl: PROGRAM IMPORTS Vector, Complex EXPORTS Complex =
BEGIN OPEN Complex;
Exponent:
PROCEDURE [x:
REAL]
RETURNS [
INTEGER] =
INLINE
BEGIN
fl: Ieee.SingleReal ← LOOPHOLE[x];
RETURN[fl.exp];
END;
AlmostEqual:
PUBLIC
PROCEDURE [a: Vec, b: Vec, mag:[-126..0] ← -20]
RETURNS [BOOLEAN] =
BEGIN
SumSqrAbs: REAL ← SqrAbs[a] + SqrAbs[b];
SqrAbsDif: REAL ← SqrAbs[Sub[a,b]];
RETURN [Exponent[SumSqrAbs]+mag+mag-1>Exponent[SqrAbsDif]];
END;
SqRt:
PUBLIC
PROCEDURE [a: Vec]
RETURNS [Vec] =
BEGIN -- Find the complex square root, using Newton's iteration
z:Vec ← (IF a.y>=0 THEN [1,1] ELSE [1,-1]);
oldz:Vec;
IF a.y = 0
THEN
BEGIN
IF a.x>0 THEN z.y𡤀 -- real square root
ELSE IF a.x<0 THEN z.x𡤀 -- pure imaginary
ELSE RETURN[[0,0]]
END;
FOR
I:
NAT
IN [0..50)
DO
oldz ← z;
z ← Vector.Mul[Add[z,Div[a,z]],0.5];
IF AlmostEqual[z,oldz,-20] THEN EXIT;
ENDLOOP;
RETURN[z];
END;
END.