<> <> <> <> <> DIRECTORY Ieee, Real, RealFns, Vector, Vector2, Complex; ComplexImpl: CEDAR PROGRAM IMPORTS RealFns, Vector, Complex EXPORTS Complex = BEGIN VEC: TYPE = Vector2.VEC; Mul: PUBLIC PROCEDURE [a: VEC, b: VEC] RETURNS [VEC] = {RETURN[[(a.x*b.x - a.y*b.y), (a.x*b.y + a.y*b.x)]]}; -- complex product Div: PUBLIC PROCEDURE [a: VEC, b: VEC] RETURNS [VEC] = BEGIN d:REAL = b.x*b.x + b.y*b.y; RETURN[[(a.x*b.x + a.y*b.y)/d, (-a.x*b.y + a.y*b.x)/d]] -- complex quotient END; FromPolar: PUBLIC PROCEDURE [r: REAL, radians: REAL] RETURNS [VEC] = {RETURN[[r*RealFns.Cos[radians], r*RealFns.Sin[radians]]]}; 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 _ Complex.SqrAbs[a] + Complex.SqrAbs[b]; sqrAbsDif: REAL _ Complex.SqrAbs[Complex.Sub[a, b]]; RETURN [Exponent[sumSqrAbs]+mag+mag-1>Exponent[sqrAbsDif]]; END; Arg: PUBLIC PROCEDURE [a: VEC] RETURNS [REAL] = {RETURN[RealFns.ArcTan[a.y, a.x]]}; Exp: PUBLIC PROCEDURE [a: VEC] RETURNS [VEC] = {RETURN[FromPolar[RealFns.Exp[a.x], a.y]]}; <> Ln: PUBLIC PROCEDURE [a: VEC] RETURNS [VEC] = {RETURN[[RealFns.Ln[Complex.Abs[a]], Arg[a]]]}; Sqr: PUBLIC PROCEDURE [a: VEC] RETURNS [VEC] = {RETURN[[(a.x*a.x - a.y*a.y), (2*a.x*a.y)]]}; 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[Complex.Add[z, Div[a, z]], 0.5]; IF AlmostEqual[z, oldz, -20] THEN EXIT; ENDLOOP; RETURN[z]; END; END.