<> <> <> 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_0 -- real square root ELSE IF a.x<0 THEN z.x_0 -- 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.