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]];
};
Init[];
END.