Ref: TYPE = Polynomial.Ref;
PolynomialRec: TYPE = Polynomial.PolynomialRec;
RealRootRec: TYPE = Polynomial.RealRootRec;
ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec;
DivideResult: TYPE = Polynomial.DivideResult;
Create:
PUBLIC
PROC [maxDegree:
NAT]
RETURNS [f:Ref] =
BEGIN
f ¬ NEW[PolynomialRec[maxDegree+1]];
FOR i:NAT IN [0..maxDegree] DO f[i] ¬ 0 ENDLOOP
END;
Degree:
PUBLIC
PROC [poly: Ref]
RETURNS [
NAT] =
BEGIN
FOR i:
NAT
DECREASING
IN [0..poly.maxDegreePlusOne)
DO
IF poly[i] # 0 THEN RETURN[i];
ENDLOOP;
RETURN[0];
END;
Eval:
PUBLIC
PROC [f: Ref, x:
REAL]
RETURNS [y:
REAL] =
BEGIN
y ¬ 0.0;
FOR i:
NAT
DECREASING
IN [0..Degree[f]]
DO
y ¬ y*x + f[i]
ENDLOOP
END;
Copy:
PUBLIC
PROC [dest, source: Ref] =
BEGIN
d: NAT ¬ Degree[source];
FOR i:NAT IN [0..d] DO dest[i] ¬ source[i] ENDLOOP;
FOR i:NAT IN (d..dest.maxDegreePlusOne) DO dest[i] ¬ 0.0 ENDLOOP;
END;
Constant:
PUBLIC
PROC [c:
ARRAY [0..0]
OF
REAL]
RETURNS [f: Ref] =
BEGIN
f ¬ Create[0];
f[0] ¬ c[0];
END;
Linear:
PUBLIC
PROC [c:
ARRAY [0..1]
OF
REAL]
RETURNS [f: Ref] =
BEGIN
f ¬ Create[1];
FOR i:NAT IN [0..1] DO f[i] ¬ c[i] ENDLOOP;
END;
Quadratic:
PUBLIC
PROC [c:
ARRAY [0..2]
OF
REAL]
RETURNS [f: Ref] =
BEGIN
f ¬ Create[2];
FOR i:NAT IN [0..2] DO f[i] ¬ c[i] ENDLOOP;
END;
Cubic:
PUBLIC
PROC [c:
ARRAY [0..3]
OF
REAL]
RETURNS [f: Ref] =
BEGIN
f ¬ Create[3];
FOR i:NAT IN [0..3] DO f[i] ¬ c[i] ENDLOOP;
END;
Quartic:
PUBLIC
PROC [c:
ARRAY [0..4]
OF
REAL]
RETURNS [f: Ref] =
BEGIN
f ¬ Create[4];
FOR i:NAT IN [0..4] DO f[i] ¬ c[i] ENDLOOP;
END;
Quintic:
PUBLIC
PROC [c:
ARRAY [0..5]
OF
REAL]
RETURNS [f: Ref] =
BEGIN
f ¬ Create[5];
FOR i:NAT IN [0..5] DO f[i] ¬ c[i] ENDLOOP;
END;
Add:
PUBLIC
PROC [dest, source: Ref] =
{FOR i:NAT IN [0..Degree[source]] DO dest[i] ¬ dest[i] + source[i] ENDLOOP};
Sum:
PUBLIC
PROC [a, b: Ref]
RETURNS [f: Ref] =
BEGIN
f ¬ Create[MAX[Degree[a],Degree[b]]];
Copy[f,a];
Add[f,b];
END;
Subtract:
PUBLIC
PROC [dest, source: Ref] =
{FOR i:NAT IN [0..Degree[source]] DO dest[i] ¬ dest[i] - source[i] ENDLOOP};
Difference:
PUBLIC
PROC [a, b: Ref]
RETURNS [f: Ref] =
BEGIN
f ¬ Create[MAX[Degree[a],Degree[b]]];
Copy[f,a];
Subtract[f,b];
END;
Product:
PUBLIC
PROC [a, b: Ref]
RETURNS [f: Ref] =
BEGIN
da:NAT ¬ Degree[a];
db:NAT ¬ Degree[b];
f ¬ Create[da+db];
FOR i:
NAT
IN [0..da]
DO
FOR j:
NAT
IN [0..db]
DO
f[i+j] ¬ f[i+j] + a[i]*b[j]
ENDLOOP
ENDLOOP
END;
DivideByLinear:
PUBLIC
PROC [f: Ref, c:
REAL]
RETURNS [value:
REAL] =
BEGIN
d: NAT ¬ Degree[f];
value ¬ 0.0;
FOR i:
NAT
DECREASING
IN [0..d]
DO
t: REAL ¬ f[i];
f[i] ¬ value;
value ¬ value*c + t
ENDLOOP;
END;
Differentiate:
PUBLIC
PROC [poly: Ref] =
BEGIN
d:NAT ¬ Degree[poly];
FOR i: NAT IN (0..d] DO poly[i-1] ¬ i*poly[i] ENDLOOP;
poly[d] ¬ 0.0;
END;
Derivative:
PUBLIC
PROC [poly: Ref]
RETURNS [dpoly: Ref] =
BEGIN
d:NAT ¬ Degree[poly];
dpoly ¬ Create[IF d>0 THEN d-1 ELSE 0];
FOR i: NAT IN (0..d] DO dpoly[i-1] ¬ i*poly[i] ENDLOOP;
END;
EvalDeriv:
PUBLIC
PROC [f: Ref, x:
REAL]
RETURNS [y:
REAL] =
BEGIN
y ¬ 0.0;
FOR i:
NAT
DECREASING
IN (0..Degree[f]]
DO
y ¬ i*f[i] + y*x
ENDLOOP
END;
Integrate:
PUBLIC
PROC [poly: Ref] =
BEGIN
d:NAT ¬ Degree[poly];
FOR i: NAT DECREASING IN (0..d+1] DO poly[i] ¬ poly[i-1]/i ENDLOOP;
poly[0] ¬ 0.0;
END;
Integral:
PUBLIC
PROC [poly: Ref]
RETURNS [f: Ref] =
BEGIN
d:NAT ¬ Degree[poly]+1;
f ¬ Create[d];
FOR i: NAT IN (0..d] DO f[i] ¬ poly[i-1]/i ENDLOOP;
END;
EvalIntegral:
PUBLIC
PROC [f: Ref, x:
REAL]
RETURNS [y:
REAL] =
BEGIN
y ¬ 0.0;
FOR i:
NAT
DECREASING
IN [0..Degree[f]]
DO
y ¬ (f[i]/(i+1) + y)*x
ENDLOOP
END;
RootBound:
PROC [poly: Ref]
RETURNS [
REAL] =
BEGIN -- finds a bound for the absolute values of the roots
d: NAT ¬ Degree[poly];
b: REAL ¬ 0;
FOR i:NAT IN [0..d) DO b ¬ b + ABS[poly[i]] ENDLOOP;
RETURN [b*1.01/ABS[poly[d]]+1];
END;
debug Bisection: SIGNAL = CODE;
RootBetween:
PROC [poly: Ref, x0, x1:
REAL]
RETURNS [x: REAL] =
BEGIN
debug xi: ARRAY [0..50) OF REAL;
debug i: NAT ← 0;
y0: REAL ¬ Eval[poly,x0];
y1: REAL ¬ Eval[poly,x1];
xx: REAL;
xx0: REAL ¬ IF y0<y1 THEN x0 ELSE x1;
xx1: REAL ¬ IF y0<y1 THEN x1 ELSE x0;
yy0: REAL ¬ MIN[y0,y1];
yy1: REAL ¬ MAX[y0,y1];
NewtonStep:
PROC [x:
REAL]
RETURNS [newx:
REAL] =
BEGIN
y: REAL ¬ Eval[poly, x];
deriv: REAL ¬ EvalDeriv[poly, x];
IF y<0 THEN {IF y>yy0 THEN {xx0 ¬ x; yy0 ¬ y}};
IF y>0 THEN {IF y<yy1 THEN {xx1 ¬ x; yy1 ¬ y}};
newx ¬ IF deriv=0 THEN x+y ELSE x-y/deriv;
END;
IF y0=0 THEN RETURN[x0];
IF y1=0 THEN RETURN[x1];
IF (y0 > 0 AND y1 > 0) OR (y0 < 0 AND y1 < 0) THEN RETURN [x1+100];
xx ¬ x ¬ (x0+x1)/2;
first try Newton's method
WHILE x
IN [x0..x1]
DO
newx:REAL ¬ NewtonStep[x];
IF x=newx THEN RETURN;
x ¬ NewtonStep[newx];
xx ¬ NewtonStep[xx];
debug xi[i] ← xx; i ← i+1;
IF xx=x THEN EXIT;
ENDLOOP;
If it oscillated or went out of range, try bisection
debug SIGNAL Bisection;
IF xx0
IN [x0..x1]
AND xx1
IN [x0..x1]
THEN
BEGIN
x0 ¬ MIN[xx0,xx1];
x1 ¬ MAX[xx0,xx1];
END;
y0 ¬ Eval[poly,x0];
y1 ¬ Eval[poly,x1];
THROUGH [0..500)
DO
y: REAL ¬ Eval[poly, x ¬ (x0+x1)/2];
IF x=x0
OR x=x1
THEN
{IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]};
IF (y > 0 AND y0 < 0) OR (y < 0 AND y0 > 0) THEN {x1 ¬ x; y1 ¬ y}
ELSE {x0 ¬ x; y0 ¬ y};
IF (y0 > 0 AND y1 > 0) OR (y0 < 0 AND y1 < 0) OR (y0=0 OR y1=0) THEN {IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]}
ENDLOOP;
ERROR;
END;
RealRoots:
PUBLIC
PROC [poly: Ref]
RETURNS [roots:
REF RealRootRec] =
BEGIN
d: NAT ¬ Degree[poly];
roots ¬ NEW[RealRootRec[d]];
roots.nRoots ¬ 0;
IF d<=1
THEN
BEGIN
IF d=1 THEN {roots.nRoots ¬ 1; roots[0] ¬ -poly[0]/poly[1]}
END
ELSE
BEGIN
extrema: REF RealRootRec ¬ RealRoots[Derivative[poly]];
bound: REAL ¬ RootBound[poly];
x: REAL;
IF extrema.nRoots>0
THEN
BEGIN
x ¬ RootBetween[poly,-bound,extrema[0]];
IF x<=extrema[0] THEN {roots.nRoots ¬ 1; roots[0] ¬ x};
FOR i:
NAT
IN [0..extrema.nRoots-1)
DO
x ¬ RootBetween[poly,extrema[i],extrema[i+1]];
IF x<=extrema[i+1] THEN {OPEN roots; roots[nRoots] ¬ x; nRoots ¬ nRoots + 1};
ENDLOOP;
x ¬ RootBetween[poly,extrema[extrema.nRoots-1],bound];
IF x<=bound THEN {OPEN roots; roots[nRoots] ¬ x; nRoots ¬ nRoots + 1}
END
ELSE
BEGIN
x ¬ RootBetween[poly,-bound,bound];
IF x<=bound THEN {OPEN roots; roots[0] ¬ x; nRoots ¬ 1}
END
END
END;
CheapRealRoots:
PUBLIC
PROC [poly: Ref]
RETURNS [roots: ShortRealRootRec] =
BEGIN
d: NAT ¬ Degree[poly];
roots.nRoots ¬ 0;
IF d<=1
THEN
BEGIN
IF d=1 THEN {roots.nRoots ¬ 1; roots.realRoot[0] ¬ -poly[0]/poly[1]}
END
ELSE
BEGIN
bound: REAL ¬ RootBound[poly];
savedCoeff: ARRAY [0..5] OF REAL;
FOR i: NAT IN [0..d] DO savedCoeff[i] ¬ poly[i] ENDLOOP;
Differentiate[poly];
BEGIN
extrema: ShortRealRootRec ¬ CheapRealRoots[poly];
x: REAL;
FOR i: NAT IN [0..d] DO poly[i] ¬ savedCoeff[i] ENDLOOP;
IF extrema.nRoots>0
THEN
BEGIN
x ¬ RootBetween[poly,-bound,extrema.realRoot[0]];
IF x<=extrema.realRoot[0] THEN {roots.nRoots ¬ 1; roots.realRoot[0] ¬ x};
FOR i:
NAT
IN [0..extrema.nRoots-1)
DO
x ¬ RootBetween[poly,extrema.realRoot[i],extrema.realRoot[i+1]];
IF x<=extrema.realRoot[i+1] THEN {roots.realRoot[roots.nRoots] ¬ x; roots.nRoots ¬ roots.nRoots + 1};
ENDLOOP;
x ¬ RootBetween[poly,extrema.realRoot[extrema.nRoots-1],bound];
IF x<=bound THEN {roots.realRoot[roots.nRoots] ¬ x; roots.nRoots ¬ roots.nRoots + 1}
END
ELSE
BEGIN
x ¬ RootBetween[poly,-bound,bound];
IF x<=bound THEN {roots.realRoot[0] ¬ x; roots.nRoots ¬ 1}
END
END
END
END;
END.
Michael Plass, June 14, 1983 3:50 pm. Fixed bug in RootBound [missing ABS on leading term]
Michael Plass, August 18, 1983 3:38 pm: Added CheapRealRoots.
--------------------------------------------------------------------------------
RealRoots: PUBLIC PROC [poly: Ref] RETURNS [roots: REF RealRootRec] =
BEGIN -- Faster, but might not always find all multiple roots.
root: REAL;
Newton: PROC [f: Ref, approx: REAL] RETURNS [BOOLEAN] =
BEGIN
root ← approx;
THROUGH [1..200] DO
fx: REAL ← Eval[f,root];
dfx: REAL ← EvalDeriv[f,root];
IF fx=0 THEN RETURN [TRUE];
IF ABS[fx]*1.0e+6 <= ABS[root*dfx] THEN RETURN [TRUE];
IF dfx = 0 THEN root ← 1-100*root
ELSE root ← root - fx/dfx
ENDLOOP;
RETURN [FALSE]
END;
d: NAT ← Degree[poly];
deflatedPoly: Ref ← Create[d];
roots ← NEW[RealRootRec[d]];
roots.nRoots ← 0;
Copy [deflatedPoly, poly];
THROUGH [1..d] WHILE Newton[deflatedPoly,0] DO
[] ← Newton[poly,root];
roots[roots.nRoots] ← root;
roots.nRoots ← roots.nRoots + 1;
[] ← DivideByLinear[deflatedPoly,root]
ENDLOOP;
END;