<> <> <> <> <> <> DIRECTORY Ieee, 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]]]}; Abs: PUBLIC PROCEDURE [a: VEC] RETURNS [REAL] = {RETURN[RealFns.SqRt[a.x*a.x+a.y*a.y]]}; -- same as Vector.Mag Exponent: PROCEDURE [x: REAL] RETURNS [INTEGER] = INLINE BEGIN RETURN[LOOPHOLE[x, Ieee.SingleReal].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.