<> <> <> <> <<>> 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; <> 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] = { <> 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] = { <> 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] = { <> 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] = { <> 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] = { <> 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] = { <> 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] = { <> 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] = { <> 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]; <> 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]; <> 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]; <> 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]; <> 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.