ComplexImpl.mesa
Copyright Ó 1985, 1989, 1992 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
Doug Wyatt, September 13, 1989 11:35:16 am PDT
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]]};
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¬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.