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 =
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.