PolynomialImpl.mesa
Michael Plass, August 19, 1983 10:17 am
DIRECTORY Polynomial;
PolynomialImpl: PROGRAM EXPORTS Polynomial =
BEGIN
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*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*y0 < 0 THEN {x1 ← x; y1 ← y}
ELSE {x0 ← x; y0 ← y};
IF y0*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;