DIRECTORY Polynomial; PolynomialImpl: CEDAR 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; RootBetween: PROC [poly: Ref, x0, x1: REAL] RETURNS [x: REAL] = BEGIN y0: REAL _ Eval[poly,x0]; y1: REAL _ Eval[poly,x1]; xx: REAL; xx0: REAL _ IF y0yy0 THEN {xx0 _ x; yy0 _ y}}; IF y>0 THEN {IF y 0 AND y1 > 0) OR (y0 < 0 AND y1 < 0) THEN RETURN [x1+100]; xx _ x _ (x0+x1)/2; WHILE x IN [x0..x1] DO newx:REAL _ NewtonStep[x]; IF x=newx THEN RETURN; x _ NewtonStep[newx]; xx _ NewtonStep[xx]; IF xx=x THEN EXIT; ENDLOOP; 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. АPolynomialImpl.mesa Copyright c 1985 by Xerox Corporation. All rights reserved. Michael Plass, July 31, 1986 3:44:22 pm PDT Tim Diebert May 21, 1985 5:52:53 pm PDT debug Bisection: SIGNAL = CODE; debug xi: ARRAY [0..50) OF REAL; debug i: NAT _ 0; first try Newton's method debug xi[i] _ xx; i _ i+1; If it oscillated or went out of range, try bisection debug SIGNAL Bisection; 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; Э╩°ШЪЬЩIcodeЪЬ ╧mЬ1ЩJЪЬЮЬЩЪ ЯЬЮЬЮЬЮЬЮЬЩ7JЪЮЩJЩЪЮЬ ЮЩJЪЬЮЬЩJЪЬЮЬЩJЪЮЬЮЬЮЬЮЬЩJЪЮЬЮЬЮЬЮЬЮЬЮЬЩ6JЪЮЬ ЮЬЩ!JЪЮЬЩJЪЮЬЩЧJЪЮЬЮЬЩJЪЮЬЩЧJЪЬЮЬЩJЩJЪЬЮЬЩJЩJЩЪЮЬЮЬЮЩ.JЩJЩJЩ JЩ&JЪЮЬЩЧJЪЮЬЩJЩЧЧЧЕЧ0z