ComplexImpl.mesa
Last edited by Michael Plass, January 4, 1984 12:00 pm
Written by Michael Plass, 29-Sep-81
DIRECTORY Ieee, Real, RealFns, Vector, Complex;
ComplexImpl: 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]]]};
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.