ComplexesImpl:
CEDAR
PROGRAM
IMPORTS IO, Convert, Rope, Real, RealFns
EXPORTS Complexes
Types and Constants
ComplexesError: PUBLIC SIGNAL [reason: ATOM ← $Unspecified] = CODE;
bitsPerWord: CARDINAL = Basics.bitsPerWord;
CARD: TYPE = LONG CARDINAL;
ROPE: TYPE = Rope.ROPE;
STREAM: TYPE = IO.STREAM;
ComplexOne: Complex ← FromPairREAL[1.0, 0.0];
ComplexZero: Complex ← FromPairREAL[0.0, 0.0];
Variables
ClassPrintName:
AC.PrintNameProc = {
RETURN["Complexes"];
};
ClassLegalFirstChar:
AC.LegalFirstCharOp = {
SELECT char
FROM
= '( => RETURN[TRUE];
ENDCASE;
RETURN[FALSE];
};
ClassRead:
AC.ReadOp = {
RETURN[Read[in]];
};
ClassFromRope:
AC.FromRopeOp = {
stream: IO.STREAM ← IO.RIS[in];
RETURN[ ClassRead[stream] ];
};
ClassToRope:
AC.ToRopeOp = {
RETURN[ ToRope[NARROW[in, Complex] ] ]
};
ClassWrite:
AC.WriteOp = {
Write[stream, NARROW[in, Complex] ]
};
ClassCharacteristic:
AC.StructureRankOp = {
RETURN[ 0 ]
};
ClassAdd:
AC.BinaryOp = {
RETURN[ Add[NARROW[firstArg, Complex], NARROW[secondArg, Complex] ] ]
};
ClassNegate:
AC.UnaryOp = {
RETURN[ Negate[ NARROW[arg, Complex] ] ];
};
ClassSubtract:
AC.BinaryOp = {
RETURN[ Subtract[NARROW[firstArg, Complex], NARROW[secondArg, Complex] ] ]
};
ClassZero:
AC.NullaryOp = {
RETURN[ ComplexZero ]
};
ClassMultiply:
AC.BinaryOp = {
RETURN[ Multiply[NARROW[firstArg, Complex], NARROW[secondArg, Complex] ] ]
};
ClassInvert:
AC.UnaryOp = {
RETURN[ Invert[NARROW[arg, Complex] ] ]
};
ClassDivide:
AC.BinaryOp = {
RETURN[ Divide[NARROW[firstArg, Complex], NARROW[secondArg, Complex] ] ]
};
ClassOne:
AC.NullaryOp = {
RETURN[ ComplexOne ]
};
ClassEqual:
AC.EqualityOp = {
RETURN[ AlmostEqual[NARROW[firstArg, Complex], NARROW[secondArg, Complex] ] ]
};
ComplexClass:
AC.StructureClass ←
NEW[
AC.StructureClassRec ← [
flavor: field,
printName: ClassPrintName,
characteristic: ClassCharacteristic,
legalFirstChar: ClassLegalFirstChar,
read: ClassRead,
fromRope: ClassFromRope,
toRope: ClassToRope,
write: ClassWrite,
add: ClassAdd,
negate: ClassNegate,
subtract: ClassSubtract,
zero: ClassZero,
multiply: ClassMultiply,
commutative: TRUE,
invert: ClassInvert,
divide: ClassDivide,
one: ClassOne,
equal: ClassEqual,
completeField: TRUE,
realField: FALSE,
realClosedField: FALSE,
algebraicallyClosedField: TRUE,
propList: NIL
] ];
Complexes:
PUBLIC
AC.Structure ←
NEW[
AC.StructureRec ← [
class: ComplexClass,
instanceData: NIL
] ];
I/O and Conversion
Read:
PUBLIC
PROC [in:
IO.
STREAM]
RETURNS [out: Complex] ~ {
puncChar, imagChar: CHAR;
var: Rope.ROPE;
real, imag: REAL;
negative: BOOL ← FALSE;
ReadVLFail: PUBLIC ERROR [subclass: ATOM ← $Unspecified] = CODE;
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[ ];
IF puncChar # '( THEN ReadVLFail[$LeftParenExpected];
[]← in.SkipWhitespace[];
real ← in.GetReal[];
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[];
IF puncChar = '- THEN negative ← TRUE ELSE IF puncChar # '+ THEN ReadVLFail[$PlusExpected];
[]← in.SkipWhitespace[];
imagChar ← in.PeekChar[];
IF imagChar # '. AND NOT imagChar IN ['0..'9] THEN imag ← 1.0 ELSE imag ← in.GetReal[];
IF negative THEN imag ← - imag;
var ← in.GetID[];
IF NOT Rope.Equal[var,"i"] AND NOT Rope.Equal[var,"I"] THEN ReadVLFail[$LetterIExpected];
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[];
IF puncChar # ') THEN ReadVLFail[$RightParenExpected];
RETURN[ FromPairREAL[real, imag] ];
};
FromRope:
PUBLIC
PROC [in: Rope.
ROPE]
RETURNS [Complex] = {
stream: IO.STREAM ← IO.RIS[in];
RETURN[ Read[ stream ] ];
};
ToRope:
PUBLIC
PROC [in: Complex]
RETURNS [out:
ROPE] = {
out ← "(";
out ← Rope.Concat[ out, RopeFromReal[in.x] ];
IF in.y < 0 THEN out ← Rope.Concat[ out, " - " ] ELSE out ← Rope.Concat[ out, " + " ];
out ← Rope.Concat[ out, RopeFromReal[ABS[in.y]] ];
out ← Rope.Concat[ out, " i )" ];
};
Write:
PUBLIC
PROC [stream:
IO.
STREAM, in: Complex] = {
IO.PutRope[ stream, ToRope[in] ]
};
FromPairREAL:
PUBLIC
PROC [realPart, imagPart:
REAL]
RETURNS [out: Complex] = {
RETURN[ NEW[ComplexRep ← [val: [realPart, imagPart] ] ] ];
};
ToPairREAL:
PUBLIC
PROC [in: Complex]
RETURNS [realPart, imagPart:
REAL] = {
RETURN[realPart: in.val.x, imagPart: in.val.y];
};
Arithmetic
Add:
PUBLIC
PROC [a, b: Complex]
RETURNS [Complex] ~ {
RETURN[NEW[ComplexRep ← [val: [a.val.x+b.val.x, a.val.y+b.val.y]] ]]
};
Negate:
PUBLIC
PROC [arg: Complex]
RETURNS [Complex] ~ {
RETURN[NEW[ComplexRep ← [val: [-arg.val.x,-arg.val.y]] ]]; -- same as vector complement
};
Subtract:
PUBLIC
PROC [a, b: Complex]
RETURNS [Complex] ~ {
RETURN[NEW[ComplexRep ← [val: [a.val.x-b.val.x, a.val.y-b.val.y]] ]]; -- same as vector difference
};
Multiply:
PUBLIC
PROC [a, b: Complex]
RETURNS [result: Complex] ~ {
RETURN[NEW[ComplexRep ← [val: [(a.val.x*b.val.x - a.val.y*b.val.y), (a.val.x*b.val.y + a.val.y*b.val.x)]] ]]; -- complex product
};
Conjugate:
PUBLIC
PROC [a: Complex]
RETURNS [Complex] ~ {
RETURN[NEW[ComplexRep ← [val: [a.val.x,-a.val.y]] ]]}; -- complex conjugate
Modulus:
PUBLIC
PROCEDURE [a: Complex]
RETURNS [
REAL] =
{RETURN[Real.SqRt[a.val.x*a.val.x+a.val.y*a.val.y]]}; -- same as Vector.Mag
ModulusSquared:
PUBLIC
PROC [a: Complex]
RETURNS [
REAL] ~ {
RETURN[a.val.x*a.val.x + a.val.y*a.val.y];
};
Invert:
PUBLIC
PROC [a: Complex]
RETURNS [result: Complex] ~ {
b: Complex ← Conjugate[a];
m: REAL ← ModulusSquared[a];
RETURN[NEW[ComplexRep ← [val: [b.val.x /m, - b.val.y /m]] ]] -- complex inverse
};
Divide:
PUBLIC
PROC [a, b: Complex]
RETURNS [result: Complex]~ {
RETURN[ Multiply[a, Invert[b]] ];
};
Exponent:
PROCEDURE [x:
REAL]
RETURNS [
INTEGER] =
INLINE
BEGIN
fl: Ieee.SingleReal ← LOOPHOLE[x];
RETURN[fl.exp];
END;
AlmostEqual:
PUBLIC
PROCEDURE [a: Complex, b: Complex, mag:[-126..0] ← -20]
RETURNS [BOOLEAN] =
BEGIN
sumSqrAbs: REAL ← ModulusSquared[a] + ModulusSquared[b];
sqrAbsDif: REAL ← ModulusSquared[Subtract[a, b]];
RETURN [Exponent[sumSqrAbs]+mag+mag-1>Exponent[sqrAbsDif]];
END;
FromPolar:
PUBLIC
PROCEDURE [r:
REAL, radians:
REAL]
RETURNS [Complex] ~
{RETURN[NEW[ComplexRep ← [val: [r*RealFns.Cos[radians], r*RealFns.Sin[radians]]]]]};
Arg:
PUBLIC
PROCEDURE [a: Complex]
RETURNS [
REAL] =
{RETURN[RealFns.ArcTan[a.val.y, a.val.x]]};
Exp:
PUBLIC
PROCEDURE [a: Complex]
RETURNS [Complex] =
{RETURN[FromPolar[RealFns.Exp[a.val.x], a.val.y]]};
complex exponential function
Ln:
PUBLIC
PROCEDURE [a: Complex]
RETURNS [Complex] =
{RETURN[NEW[ComplexRep ← [val: [RealFns.Ln[Modulus[a]], Arg[a]]]]]};
Sqr:
PUBLIC
PROCEDURE [a: Complex]
RETURNS [Complex] =
{RETURN[NEW[ComplexRep ← [val: [(a.val.x*a.val.x - a.val.y*a.val.y), (2*a.val.x*a.val.y)]]]]};
SqRt:
PUBLIC
PROCEDURE [a: Complex]
RETURNS [Complex] =
BEGIN -- Find the complex square root, using Newton's iteration
z: Complex ← (IF a.val.y>=0 THEN FromPairREAL[1, 1] ELSE FromPairREAL[1, -1]);
oldz:Complex;
IF a.val.y = 0
THEN
BEGIN
IF a.val.x>0 THEN z.y𡤀 -- real square root
ELSE IF a.val.x<0 THEN z.x𡤀 -- pure imaginary
ELSE RETURN[FromPairREAL[0, 0]]
END;
FOR
I:
NAT
IN [0..50)
DO
oldz ← z;
z ← Add[z, Divide[a, z]];
z ← NEW[ComplexRep ← [val: [z.val.x / 0.5, z.val.y / 0.5 ] ] ];
IF AlmostEqual[z, oldz, -20] THEN EXIT;
ENDLOOP;
RETURN[z];
END;