ComplexImpl.mesa
Copyright © 1985 by Xerox Corporation. All rights reserved.
Last edited by Michael Plass, January 4, 1984 12:00 pm
Written by Michael Plass, 29-Sep-81
Tim Diebert May 20, 1985 2:07:16 pm PDT
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]]};
complex exponential function
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𡤀 -- 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[Complex.Add[z, Div[a, z]], 0.5];
IF AlmostEqual[z, oldz, -20] THEN EXIT;
ENDLOOP;
RETURN[z];
END;
END.