DIRECTORY Rope, IO, Basics, Ieee, Real, RealFns, Vector, Vector2, Atom, Convert, AlgebraClasses, Points, Complexes; ComplexesImpl: CEDAR PROGRAM IMPORTS IO, Convert, Rope, Real, RealFns EXPORTS Complexes = BEGIN OPEN Complexes, Convert, AC: AlgebraClasses, PTS: Points; 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]; 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 ] ]; 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]; }; 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]]}; 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_0 -- real square root ELSE IF a.val.x<0 THEN z.x_0 -- 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; END. ¨ComplexesImpl.mesa Last Edited by: Arnon, July 19, 1985 2:54:46 pm PDT Types and Constants Variables I/O and Conversion Arithmetic complex exponential function Κ Ψ˜Jšœ™J™3J™šΟk ˜ Jšœ˜Icodešœ˜Kšœ˜Jšœ˜K˜Kšœ˜Kšœ˜Kšœ ˜ K˜K˜J˜J˜Jšœ ˜ —J˜head2šœœ˜Jšœœ˜(Jšœ ˜J˜—Jš œœœœœ ˜AheadšΟn™Jš œœœ œœ˜CKšœ œ˜+Kšœœœœ˜Kšœœœ˜Kšœœœœ˜Kšœ-˜-Kšœ.˜.—šž ™ šžœœ˜$Jšœ˜J˜—šžœœ˜,šœ˜Kšœœœ˜Kšœ˜—Kšœœ˜J˜—šž œœ ˜Jšœ ˜J˜—šž œœ˜ Kš œœœœœ˜Kšœ˜K˜—šž œœ ˜Jšœ œ˜&J˜—šž œœ ˜Jšœœ˜#Jšœ˜—šžœœ˜+Jšœ˜ Jšœ˜—šžœœ ˜Jšœœœ˜EJšœ˜—šž œœ ˜Jšœ œ˜)Jšœ˜—šž œœ ˜Jšœ œœ˜JJšœ˜—šž œœ˜Jšœ˜Jšœ˜—šž œœ ˜Jšœ œœ˜JJšœ˜—šž œœ ˜Jšœ œ˜'Jšœ˜—šž œœ ˜Jšœ œœ˜HJšœ˜—šžœœ˜Jšœ˜Jšœ˜—šž œœ˜Jšœœœ˜MJšœ˜—J˜šž œœœœ˜>J˜J˜Jšœ$˜$J˜Jšœ$˜$Jšœ˜Jšœ˜Jšœ˜Jšœ˜J˜Kšœ˜Kšœ˜Kšœ˜K˜K˜Kšœ˜Kšœ œ˜Kšœ˜Kšœ˜K˜K˜K˜K˜Kšœœ˜Kšœ œ˜Kšœœ˜Kšœœ˜K˜Jšœ ˜ K˜—K˜š ž œœœ œœ˜8Jšœ˜Jšœ˜K˜——šœ™š žœœœœœœ˜Jšœ˜Jšœœ˜Jšœœ3Ÿ˜OK˜J˜—šžœœœœ˜@Jšœ˜!K˜—š žœ œœœœ˜8Jš˜Jšœœ˜"Jšœ ˜Jšœ˜J˜—šž œœ œ.˜KJšœœ˜Jš˜Jšœ œ)˜8Jšœ œžœ˜1Jšœ5˜;Jšœ˜—š ž œœ œœ œœ ˜HJšœœœI˜TJ˜—š žœœ œœœ˜3Jšœœ$˜+J˜—šžœœ œœ ˜6Jšœœ,˜3Jšœ™J˜—šžœœ œœ ˜5Jšœœœ9˜DJ˜—šžœœ œœ ˜6JšœœœS˜^J˜—šžœœ œœ ˜7JšœŸ9˜?Jšœœ œœ˜NJšœ ˜ šœ œ˜Jš˜Jšœ œŸ˜+Jšœœ œŸ˜.Jšœœ˜Jšœ˜—š œžœœœ ˜J˜ Jšœ žœ˜Jšœœ8˜?Jšœœœ˜'Jšœ˜—Jšœ˜ Jšœ˜—J˜˜J˜J˜J˜——˜J˜—Jšœ˜J˜—…— 'Š