File: RootsImpl.mesa
Copyright © 1984 by Xerox Corporation. All rights reserved.
Last edited by Bier on August 1, 1985 1:42:01 am PDT
Contents: The quadratic formula, cubic formula and quartic formula as translated from LISP from
"LISP" by Patrick Winston and Berthold Horn. All functions find complex roots.
DIRECTORY
Complex,
FS,
IO,
Polynomial,
RealFns,
Roots,
Rope;
RootsImpl:
PROGRAM
IMPORTS Complex, FS, IO, Polynomial, RealFns
EXPORTS Roots =
BEGIN
Vec: TYPE = Complex.VEC;
RootArray: TYPE = Roots.RootArray;
LinRoots: TYPE = Roots.LinRoots;
CubeRoots: TYPE = Roots.CubeRoots;
QuarticRoots: TYPE = Roots.QuarticRoots;
testCount: NAT = 17;
shortTestCount: NAT = 5;
TestVector: TYPE = ARRAY [1..testCount] OF REAL;
ShortTestVector: TYPE = ARRAY [1..shortTestCount] OF REAL;
GLOBALS
Inconsistent: PUBLIC SIGNAL = CODE;
root3over2: REAL ← RealFns.SqRt[3.0]/2.0;
coeffs: TestVector ← [-25.0, -17.0, -12.0, -5.6, -2.2, -1.0, -0.5, -0.001, 0.0, 0.001, 0.5, 1.0, 2.2, 5.6, 12.0, 17.0, 25.0];
shortCoeffs: ShortTestVector ← [-17.0, -1.0, 0.0, 1.0, 17.0];
ComplexLinearFormula:
PUBLIC
PROC [a, b:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
The solution to the equation "ax + b = 0".
IF a = 0
THEN {
IF b = 0 THEN {roots[1] ← [0.0, 0.0]; rootCount ← 0; RETURN}
ELSE roots[1] ← [0.0, 0.0]; rootCount ← 0; RETURN} -- inconsistent case
ELSE {
IF b = 0 THEN {roots[1] ← [0.0, 0.0]; rootCount ← 1; RETURN}
ELSE {roots[1] ← [-b/a, 0.0]; rootCount ← 1; RETURN};
};
};
ComplexQuadratic:
PUBLIC
PROC [a, b, c:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
The solution to the equation "ax2+bx+c=0". If a=0 this is just a linear equation. If
c = 0, one root is zero and we solve a linear equation. Otherwise, we use the
quadratic formula in either the form [-b+-(b2-4ac)1/2]/2a or (2c)/[-b-+(b2-4ac)1/2]
depending on ...
IF a = 0
THEN {
linRoots: RootArray;
[linRoots, rootCount] ← ComplexLinearFormula[b, c];
RETURN[linRoots, rootCount];
};
IF c = 0
THEN {
linRoots: RootArray;
roots[1] ← [0.0, 0.0];
[linRoots, rootCount] ← ComplexLinearFormula[a, b];
roots[2] ← linRoots[1];
rootCount ← rootCount + 1;
RETURN};
RETURN QuadraticAux[a, b, c, (b*b)-(4.0*a*c)];
}; -- end of ComplexQuadraticFormula
QuadraticAux:
PRIVATE
PROC [a, b, c, discriminant:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
IF discriminant < 0
THEN {
twoA: REAL ← 2.0*a;
RETURN QuadraticConjugate[[-b/twoA, RealFns.SqRt[-discriminant]/twoA]];
};
IF discriminant = 0
THEN {
twoA: REAL ← 2.0*a;
roots[1] ← roots[2] ← [-b/twoA, 0.0];
rootCount ← 2;
RETURN;
};
IF b < 0 THEN RETURN QuadraticReal[a, RealFns.SqRt[discriminant] - b, c]
ELSE RETURN QuadraticReal[a, -RealFns.SqRt[discriminant] - b, c];
}; -- end of QuadraticAux
QuadraticReal:
PRIVATE
PROC [a, term, c:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
roots[1] ← [term/(2.0*a), 0.0];
roots[2] ← [(2.0*c)/term, 0.0];
rootCount ← 2;
RETURN;
}; -- end of QuadraticReal
QuadraticConjugate:
PRIVATE
PROC [z: Vec]
RETURNS [roots: RootArray, rootCount:
NAT] = {
roots[1] ← z;
roots[2] ← [z.x, -z.y];
rootCount ← 2;
}; -- end of QuadraticConjugate
ComplexCubic:
PUBLIC
PROC [a, b, c, d:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
Finds the roots of equations of the form "ax^3+bx^2+cx+d".
bb, bbb, cubeTerm: REAL;
quadRoots: RootArray;
quadCount: NAT;
IF a = 0 THEN RETURN ComplexQuadratic[b, c, d];
bb ← b*b;
bbb ← bb*b;
cubeTerm ← (bb - 3.0*a*c);
cubeTerm ← cubeTerm*cubeTerm*cubeTerm;
[quadRoots, quadCount] ← ComplexQuadratic[1.0, 2.0*bbb + 9.0*a*(3.0*a*d-b*c), cubeTerm];
RETURN CubicAux[a,
b,
quadRoots,
quadCount];
}; -- end of ComplexCubic
CubicAux:
PRIVATE
PROC [a, b:
REAL, quadRoots: RootArray, quadCount:
NAT]
RETURNS [roots: RootArray, rootCount:
NAT] = {
IF quadRoots[1].y = 0
THEN
-- quadRoots real
RETURN CubicConjugate[a, b, CubeRoot[quadRoots[1].x], CubeRoot[quadRoots[2].x]]
ELSE RETURN CubicReal[a, b, RealFns.SqRt[(quadRoots[1].x*quadRoots[1].x) +
(quadRoots[1].y*quadRoots[1].y)],
RealFns.ArcTan[quadRoots[1].y, quadRoots[1].x]];
}; -- end of CubicAux
CubeRoot:
PUBLIC
PROCEDURE [z:
REAL]
RETURNS [root:
REAL] = {
IF z < 0
THEN {
posZ: REAL ← -z;
root ← RealFns.Exp[RealFns.Ln[posZ]/3.0];
root ← -(2.0*root + posZ/(root*root))/3.0;
}
ELSE {
root ← RealFns.Exp[RealFns.Ln[z]/3.0];
root ← (2.0*root + z/(root*root))/3.0;
};
};
CubicConjugate:
PRIVATE
PROC [a, b, r, s:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
r and s are cube roots
threeA: REAL ← 3.0*a;
RETURN CubicConjugateAux[(r+s-b)/threeA, (-(r+s)/2.0 - b)/threeA, (r-s)*root3over2/threeA];
}; -- end of CubicConjugate
CubicConjugateAux:
PRIVATE
PROC [realRoot, real, imaginary:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
roots[1] ← [realRoot, 0.0];
roots[2] ← [real, imaginary];
roots[3] ← [real, -imaginary];
rootCount ← 3;
}; -- end of CubicConjugateAux
CubicReal:
PRIVATE
PROC [a, b, rho, theta:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
RETURN CubicRealAux[a, b, 2.0*CubeRoot[rho], RealFns.Cos[theta/3.0]*-0.5,
RealFns.Sin[theta/3.0]*-root3over2];
};
CubicRealAux:
PROCEDURE [a, b, rd, cd, sd:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
threeA: REAL ← 3.0*a;
roots[1] ← [((-2.0*rd*cd) - b)/threeA, 0.0];
roots[2] ← [(rd*(cd+sd) - b)/threeA, 0.0];
roots[3] ← [(rd*(cd-sd) - b)/threeA, 0.0];
rootCount ← 3;
};
ComplexQuartic:
PROCEDURE [a, b, c, d, e:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
cubicCount: NAT;
cubicRoots: RootArray;
IF a=0 THEN RETURN ComplexCubic[b, c, d, e];
IF e=0
THEN {
[cubicRoots, cubicCount] ← ComplexCubic[a, b, c, d];
roots[1] ← [0.0, 0.0];
FOR i:
NAT
IN [1..cubicCount]
DO
roots[i+1] ← cubicRoots[i];
ENDLOOP;
rootCount ← cubicCount+1;
RETURN;
};
[cubicRoots, cubicCount] ← ComplexCubic[1.0, -c, b*d-4.0*a*e, 4.0*a*c*e-(a*d*d+b*b*e)];
RETURN QuarticAux[a, b, c, d, e, cubicRoots[1].x];
}; -- end of ComplexQuartic
QuarticAux:
PROCEDURE [a, b, c, d, e, s:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
RETURN QuarticSplit[a, b, RealFns.SqRt[ABS[b*b - 4.0*a*(c-s)]], s,
RealFns.SqRt[ABS[s*s-4.0*a*e]], b*s - 2.0*a*d];
}; -- end of QuarticAux
QuarticSplit:
PROCEDURE [a, b, r1, s, r2, bsMinus2ad:
REAL]
RETURNS [roots: RootArray, rootCount:
NAT] = {
quadRoots: RootArray;
quadCount: NAT;
IF (r1*r2*bsMinus2ad) < 0
THEN {
[roots, rootCount] ← ComplexQuadratic[2.0*a, b - r1, s+r2];
[quadRoots, quadCount] ← ComplexQuadratic[2.0*a, b+r1, s-r2];
roots[rootCount+1] ← quadRoots[1];
roots[rootCount+2] ← quadRoots[2];
rootCount ← rootCount + quadCount;
RETURN;
}
ELSE {
[roots, rootCount] ← ComplexQuadratic[2.0*a, b - r1, s-r2];
[quadRoots, quadCount] ← ComplexQuadratic[2.0*a, b+r1, s+r2];
roots[rootCount+1] ← quadRoots[1];
roots[rootCount+2] ← quadRoots[2];
rootCount ← rootCount + quadCount;
RETURN;
};
}; -- end of QuarticSplit
RealQuartic:
PUBLIC
PROC [a, b, c, d, e:
REAL]
RETURNS [roots:
ARRAY[1..4]
OF
REAL, rootCount:
NAT] = {
complexRoots: RootArray;
complexCount: NAT;
recipMaxMag: REAL;
recipMaxMag ← 1.0/MaxFive[a, b, c, d, e];
[complexRoots, complexCount] ← ComplexQuartic[a*recipMaxMag, b*recipMaxMag, c*recipMaxMag, d*recipMaxMag, e*recipMaxMag];
rootCount ← 0;
FOR i:
NAT
IN [1..complexCount]
DO
IF complexRoots[i].y = 0.0
THEN {
rootCount ← rootCount + 1;
roots[rootCount] ← complexRoots[i].x;
};
ENDLOOP;
}; -- end of RealQuartic
MaxFive:
PROCEDURE [a, b, c, d, e:
REAL]
RETURNS [max:
REAL] = {
aMag, bMag, cMag, dMag, eMag: REAL;
aMag ← ABS[a]; bMag ← ABS[b]; cMag ← ABS[c]; dMag ← ABS[d]; eMag ← ABS[e];
max ← aMag;
IF bMag>max THEN max ← bMag;
IF cMag>max THEN max ← cMag;
IF dMag>max THEN max ← dMag;
IF eMag>max THEN max ← eMag;
};
EvalQuartic:
PROCEDURE [r: Vec, a, b, c, d, e:
REAL]
RETURNS [result: Vec] = {
Evaluate the quartic polynomial ar4 + br3 + cr2 + dr + e using Horner's rule:
(((ar+b)r+c)r+d)r+e.
OPEN Complex;
t: Vec;
t ← ComplexAdd[ComplexMul[[a, 0.0], r], [b, 0.0]]; -- t ← ar+b;
t ← ComplexAdd[ComplexMul[t, r], [c, 0.0]]; -- t ← (ar+b)r+c;
t ← ComplexAdd[ComplexMul[t, r], [d, 0.0]]; -- t ← ((ar+b)r+c)r+d;
t ← ComplexAdd[ComplexMul[t, r], [e, 0.0]]; -- t ← (((ar+b)r+c)r+d)r+e;
result ← t;
};
EvalQuarticAtReal:
PROCEDURE [r:
REAL, a, b, c, d, e:
REAL]
RETURNS [result:
REAL] = {
Evaluate the quartic polynomial ar4 + br3 + cr2 + dr + e using Horner's rule:
(((ar+b)r+c)r+d)r+e.
OPEN Complex;
t: REAL;
t ← (((a*r+b)*r+c)*r+d)*r+e;
result ← t;
};
EvalCubic:
PROCEDURE [r: Vec, a, b, c, d:
REAL]
RETURNS [result: Vec] = {
Evaluate the quartic polynomial ar3 + br2 + cr + d using Horner's rule:
((ar+b)r+c)r+d.
OPEN Complex;
t: Vec;
t ← ComplexAdd[ComplexMul[[a, 0.0], r], [b, 0.0]]; -- t ← ar+b;
t ← ComplexAdd[ComplexMul[t, r], [c, 0.0]]; -- t ← (ar+b)r+c;
t ← ComplexAdd[ComplexMul[t, r], [d, 0.0]]; -- t ← ((ar+b)r+c)r+d;
result ← t;
};
EvalQuad:
PROCEDURE [r: Vec, a, b, c:
REAL]
RETURNS [result: Vec] = {
Evaluate the quartic polynomial ar2 + br + c using Horner's rule:
(ar+b)r+c.
OPEN Complex;
t: Vec;
t ← ComplexAdd[ComplexMul[[a, 0.0], r], [b, 0.0]]; -- t ← ar+b;
t ← ComplexAdd[ComplexMul[t, r], [c, 0.0]]; -- t ← (ar+b)r+c;
result ← t;
};
TestQuad:
PUBLIC
PROC [filename: Rope.
ROPE] = {
file: IO.STREAM;
roots: RootArray;
rootCount: NAT;
eval: Vec;
file ← FS.StreamOpen[filename, $create];
Letting coefficients a, and b take on any of the values in the test vector, we use our quad finder to find all of the complex roots of each polynomial. For each root, we evaluate the polynomial at that root. The result should come out zero. For now, write out a file entry for each root
FOR i:
NAT
IN [1..testCount]
DO
FOR j:
NAT
IN [1..testCount]
DO
[roots, rootCount] ← ComplexQuadratic[coeffs[i], coeffs[j], 1.0];
file.PutF["poly: [%g, %g, %g] ", [real[coeffs[i]]], [real[coeffs[j]]], [real[1.0]]];
file.PutF["roots: [(%g, %g), (%g, %g)]\n", [real[roots[1].x]], [real[roots[1].y]],
[real[roots[2].x]], [real[roots[2].y]]];
FOR k:
NAT
IN [1..rootCount]
DO
eval ← EvalQuad[roots[k], coeffs[i], coeffs[j], 1.0];
file.PutF[" root %g evals to (%g, %g).\n", [integer[k]], [real[eval.x]], [real[eval.y]]];
ENDLOOP;
ENDLOOP;
ENDLOOP;
file.Close[];
};
TestCubic:
PUBLIC
PROC [filename: Rope.
ROPE] = {
file: IO.STREAM;
roots: RootArray;
rootCount: NAT;
eval: Vec;
file ← FS.StreamOpen[filename, $create];
Letting coefficients a, b, and c take on any of the values in the test vector, we use our quad finder to find all of the complex roots of each polynomial. For each root, we evaluate the polynomial at that root. The result should come out zero. For now, write out a file entry for each root.
FOR i:
NAT
IN [1..shortTestCount]
DO
FOR j:
NAT
IN [1..shortTestCount]
DO
FOR l:
NAT
IN [1..shortTestCount]
DO
[roots, rootCount] ← ComplexCubic[shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], 1.0];
file.PutF["poly: [%g, %g, %g, %g] ", [real[shortCoeffs[i]]], [real[shortCoeffs[j]]], [real[shortCoeffs[l]]], [real[1.0]]];
file.PutF["roots: ["];
file.PutF["[%g, %g]", [real[roots[1].x]], [real[roots[1].y]]];
FOR k:
NAT
IN [2..rootCount]
DO
file.PutF[", [%g, %g]", [real[roots[k].x]], [real[roots[k].y]]];
ENDLOOP;
file.PutF["]\n"];
FOR k:
NAT
IN [1..rootCount]
DO
eval ← EvalCubic[roots[k], shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], 1.0];
file.PutF[" root %g evals to (%g, %g).\n", [integer[k]], [real[eval.x]], [real[eval.y]]];
ENDLOOP;
ENDLOOP;
ENDLOOP;
ENDLOOP;
file.Close[];
};
TestQuartic:
PUBLIC
PROC [filename: Rope.
ROPE] = {
file: IO.STREAM;
roots: RootArray;
rootCount: NAT;
eval: Vec;
file ← FS.StreamOpen[filename, $create];
Letting coefficients a, b, and c take on any of the values in the test vector, we use our quad finder to find all of the complex roots of each polynomial. For each root, we evaluate the polynomial at that root. The result should come out zero. For now, write out a file entry for each root.
FOR i:
NAT
IN [1..shortTestCount]
DO
FOR j:
NAT
IN [1..shortTestCount]
DO
FOR l:
NAT
IN [1..shortTestCount]
DO
FOR m:
NAT
IN [1..shortTestCount]
DO
[roots, rootCount] ← ComplexQuartic[shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
file.PutF["poly: [%g, %g, %g, %g, %g] ", [real[shortCoeffs[i]]], [real[shortCoeffs[j]]], [real[shortCoeffs[l]]], [real[shortCoeffs[m]]], [real[1.0]]];
file.PutF["roots: ["];
file.PutF["[%g, %g]", [real[roots[1].x]], [real[roots[1].y]]];
FOR k:
NAT
IN [2..rootCount]
DO
file.PutF[", [%g, %g]", [real[roots[k].x]], [real[roots[k].y]]];
ENDLOOP;
file.PutF["]\n"];
FOR k:
NAT
IN [1..rootCount]
DO
eval ← EvalQuartic[roots[k], shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
file.PutF[" root %g evals to (%g, %g).\n", [integer[k]], [real[eval.x]], [real[eval.y]]];
ENDLOOP;
ENDLOOP;
ENDLOOP;
ENDLOOP;
ENDLOOP;
file.Close[];
};
globalPoly: Polynomial.Ref;
CompareQuartic:
PUBLIC
PROC [filename: Rope.
ROPE] = {
file: IO.STREAM;
roots: ARRAY[1..4] OF REAL;
rootCount: NAT;
realRoots: REF Polynomial.RealRootRec;
problem: BOOL;
eval: REAL;
file ← FS.StreamOpen[filename, $create];
Letting coefficients a, b, and c take on any of the values in the test vector, we use our two quad finders to find all of the real roots of each polynomial. If the finders disagree on the number of roots, we print out the complete situation. For each root, we evaluate the polynomial at that root. If either finder is off by more than e-3, we print out the situation.
FOR i:
NAT
IN [1..shortTestCount]
DO
FOR j:
NAT
IN [1..shortTestCount]
DO
FOR l:
NAT
IN [1..shortTestCount]
DO
FOR m:
NAT
IN [1..shortTestCount]
DO
[roots, rootCount] ← RealQuartic[shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
globalPoly[4] ← shortCoeffs[i];
globalPoly[3] ← shortCoeffs[j];
globalPoly[2] ← shortCoeffs[l];
globalPoly[1] ← shortCoeffs[m];
globalPoly[0] ← 1.0;
realRoots ← Polynomial.RealRoots[globalPoly];
IF realRoots.nRoots # rootCount
THEN {
file.PutF["Different number of roots ** ** ** **\n"];
file.PutF["poly: [%g, %g, %g, %g, %g] ", [real[shortCoeffs[i]]], [real[shortCoeffs[j]]], [real[shortCoeffs[l]]], [real[shortCoeffs[m]]], [real[1.0]]];
file.PutF["Closed form:\n"];
file.PutF["roots: (%g)[", [integer[rootCount]] ];
IF rootCount > 0 THEN {
file.PutF["%g", [real[roots[1]]] ];
FOR k:
NAT
IN [2..rootCount]
DO
file.PutF[", %g", [real[roots[k]]] ];
ENDLOOP;
};
file.PutF["]\n"];
FOR k:
NAT
IN [1..rootCount]
DO
eval ← EvalQuarticAtReal[roots[k], shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
file.PutF[" root %g evals to %g.\n", [integer[k]], [real[eval]] ];
ENDLOOP;
file.PutF["Iterative form:\n"];
file.PutF["roots: (%g)[", [integer[realRoots.nRoots]] ];
IF realRoots.nRoots > 0 THEN {
file.PutF["%g", [real[realRoots[0]]] ];
FOR k:
NAT
IN [1..realRoots.nRoots-1]
DO
file.PutF[", %g", [real[realRoots[k]]] ];
ENDLOOP;
};
file.PutF["]\n"];
FOR k:
NAT
IN [0..realRoots.nRoots)
DO
eval ← EvalQuarticAtReal[realRoots[k], shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
file.PutF[" root %g evals to %g.\n", [integer[k]], [real[eval]] ];
ENDLOOP;
}
ELSE {
problem ← FALSE;
FOR k:
NAT
IN [1..rootCount]
DO
eval ← EvalQuarticAtReal[roots[k], shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
IF eval > 1.0e-3 THEN problem ← TRUE;
ENDLOOP;
FOR k:
NAT
IN [0..realRoots.nRoots)
DO
eval ← EvalQuarticAtReal[realRoots[k], shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
IF eval > 1.0e-3 THEN problem ← TRUE;
ENDLOOP;
IF problem
THEN {
file.PutF["EVAL TOO LARGE\n"];
file.PutF["poly: [%g, %g, %g, %g, %g] ", [real[shortCoeffs[i]]], [real[shortCoeffs[j]]], [real[shortCoeffs[l]]], [real[shortCoeffs[m]]], [real[1.0]]];
file.PutF["Closed form:\n"];
file.PutF["roots: (%g)[", [integer[rootCount]] ];
file.PutF["%g", [real[roots[1]]] ];
FOR k:
NAT
IN [2..rootCount]
DO
file.PutF[", %g", [real[roots[k]]] ];
ENDLOOP;
file.PutF["]\n"];
FOR k:
NAT
IN [1..rootCount]
DO
eval ← EvalQuarticAtReal[roots[k], shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
file.PutF[" root %g evals to %g.\n", [integer[k]], [real[eval]] ];
ENDLOOP;
file.PutF["Iterative form:\n"];
file.PutF["roots: (%g)[", [integer[realRoots.nRoots]] ];
file.PutF["%g", [real[realRoots[0]]] ];
FOR k:
NAT
IN [1..realRoots.nRoots-1]
DO
file.PutF[", %g", [real[realRoots[k]]] ];
ENDLOOP;
file.PutF["]\n"];
FOR k:
NAT
IN [0..realRoots.nRoots)
DO
eval ← EvalQuarticAtReal[realRoots[k], shortCoeffs[i], shortCoeffs[j], shortCoeffs[l], shortCoeffs[m], 1.0];
file.PutF[" root %g evals to %g.\n", [integer[k]], [real[eval]] ];
ENDLOOP;
};
};
ENDLOOP;
ENDLOOP;
ENDLOOP;
ENDLOOP;
file.Close[];
};
ComplexMul:
PROC [r, s: Vec]
RETURNS [z: Vec] = {
z ← Complex.Mul[r, s];
};
ComplexAdd:
PROC [r, s: Vec]
RETURNS [z: Vec] = {
z ← Complex.Add[r, s];
};
Init:
PROC = {
globalPoly ← Polynomial.Quartic[[0,0,0,0,0]];
};
END.