<> <> <> <> DIRECTORY SVMasterObject, SVCastRays, SVRay, SVMasterObjectTypes, SVObjectCast, Polynomial, RealFns, SV2d, SV3d, SVSceneTypes, Roots; SVObjectCastImplC: CEDAR PROGRAM IMPORTS SVCastRays, SVRay, Polynomial, RealFns EXPORTS SVObjectCast = BEGIN Classification: TYPE = SVSceneTypes.Classification; Point2d: TYPE = SV2d.Point2d; Point3d: TYPE = SV3d.Point3d; Primitive: TYPE = SVSceneTypes.Primitive; Ray: TYPE = SVSceneTypes.Ray; Surface: TYPE = SVSceneTypes.Surface; TorusRec: TYPE = SVMasterObjectTypes.TorusRec; Vector3d: TYPE = SV3d.Vector3d; <> ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec; Ref: TYPE = Polynomial.Ref; globalPoly: Polynomial.Ref; newWay: BOOL _ TRUE; <> PositiveRoots: PROCEDURE [a4, a3, a2, a1, a0: REAL] RETURNS [rootArray: Polynomial.ShortRealRootRec] = { <> i: NAT _ 1; globalPoly[4] _ a4; -- optional globalPoly[3] _ a3; -- optional globalPoly[2] _ a2; -- optional globalPoly[1] _ a1; -- optional globalPoly[0] _ a0; -- optional rootArray _ CheapRealRoots[globalPoly]; -- optional <<[realRootsOneOrigin, rootCount] _ Roots.RealQuartic[a4, a3, a2, a1, a0];>> <> <> IF rootArray.nRoots = 0 THEN RETURN; WHILE i <= rootArray.nRoots DO IF rootArray.realRoot[i-1] <= 0 THEN { FOR j: NAT IN [i..rootArray.nRoots-1] DO rootArray.realRoot[j-1] _ rootArray.realRoot[j]; ENDLOOP; rootArray.nRoots _ rootArray.nRoots - 1; } ELSE {i _ i + 1}; ENDLOOP; }; ToroidCast: PUBLIC PROC [localRay: Ray, prim: Primitive, torusRec: REF ANY, surface: Surface] RETURNS [class: Classification] = { IF newWay THEN RETURN[NewToroidCast[localRay, prim, torusRec, surface]] ELSE RETURN[OldToroidCast[localRay, prim, torusRec, surface]]; }; OldToroidCast: PUBLIC PROC [localRay: Ray, prim: Primitive, torusRec: REF ANY, surface: Surface] RETURNS [class: Classification] = { torusData: TorusRec; realRoots: Polynomial.ShortRealRootRec; -- optional <> <> <> B,C: REAL; n,N,q,Q,l,L: REAL; BB, CC: REAL; D,E,F: REAL; nn, NN, qq, QQ, ll, LL, NNNN, QQQQ, LLLL, nN, qQ, lL: REAL; a4,a3,a2,a1,a0: REAL; p: Point3d; d: Vector3d; class _ SVCastRays.GetClassFromPool[]; torusData _ NARROW[torusRec]; [p, d] _ SVRay.GetLocalRay[localRay]; B _ torusData.bigRadius; C _ torusData.sectionRadius; n _ p[1]; N _ d[1]; q _ p[2]; Q _ d[2]; l _ p[3]; L _ d[3]; BB _ B*B; CC _ C*C; D _ BB + CC; E _ BB - CC; F _ BB*CC; nn _ n*n; NN _ N*N; qq _ q*q; QQ _ Q*Q; ll _ l*l; LL _ L*L; NNNN _ NN*NN; QQQQ _ QQ*QQ; LLLL _ LL*LL; nN _ n*N; qQ _ q*Q; lL _ l*L; a4 _ NNNN + QQQQ + LLLL + 2*(NN*(QQ + LL) + QQ*LL); a3 _ 4*(nN*NN + qQ*QQ + lL*LL + nN*(QQ + LL) + NN*(qQ + lL) + lL*QQ + qQ*LL); a2 _ 2*(3*(nn*NN + qq*QQ + ll*LL) + NN*(qq + ll) + nn*(QQ+LL) + 4*nN*(qQ+lL) - D*NN + E*QQ - D*LL + ll*QQ + LL*qq + 4*qQ*lL); a1 _ 4*(nn*nN + qq*qQ + ll*lL + nN*(qq + ll) + (qQ + lL)*nn - nN*D + qQ*E - lL*D + qQ*ll + lL*qq); a0 _ nn*nn + qq*qq + ll*ll + BB*BB + CC*CC + 2*(nn*(qq + ll) - D*nn + E*qq - D*ll + qq*ll - F); realRoots _ PositiveRoots[1.0, a3/a4, a2/a4, a1/a4, a0/a4]; <<>> SELECT realRoots.nRoots FROM <> 0 => {-- complete miss. class.count _ 0; class.classifs[1] _ FALSE; }; 1 => {-- a skim. Count as a miss. class.count _ 0; class.classifs[1] _ FALSE; }; 2 => {-- cuts one cross section of the doughnut. Sort roots by size. x, y, z, r, OneMinusBoverR: REAL; class.count _ 2; IF realRoots.realRoot[0] < realRoots.realRoot[1] THEN {class.params[1] _ realRoots.realRoot[0]; class.params[2] _ realRoots.realRoot[1]} ELSE {class.params[2] _ realRoots.realRoot[0]; class.params[1] _ realRoots.realRoot[1]}; class.surfaces[1] _ surface; class.surfaces[2] _ surface; class.classifs[1] _ FALSE; class.classifs[2] _ TRUE; class.classifs[3] _ FALSE; class.primitives[1] _ class.primitives[2] _ prim; <> x _ n + class.params[1]*N; y _ q + class.params[1]*Q; z _ l + class.params[1]*L; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; <> x _ n + class.params[2]*N; y _ q + class.params[2]*Q; z _ l + class.params[2]*L; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[2] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; }; 3 => {-- cuts one cross section. Skims the other. For now, count as a miss. class.count _ 0; class.classifs[1] _ FALSE; }; 4 => {-- cuts both cross sections. Sort the roots. Pair them in order. x, y, z, r, OneMinusBoverR: REAL; <> SortHitRec[realRoots]; class.count _ 4; class.params[1] _ realRoots.realRoot[0]; class.params[2] _ realRoots.realRoot[1]; class.params[3] _ realRoots.realRoot[2]; class.params[4] _ realRoots.realRoot[3]; class.surfaces[1] _ class.surfaces[2] _ class.surfaces[3] _ class.surfaces[4] _ surface; class.classifs[1] _ FALSE; class.classifs[2] _ TRUE; class.classifs[3] _ FALSE; class.classifs[4] _ TRUE; class.classifs[5] _ FALSE; class.primitives[1] _ class.primitives[2] _ class.primitives[3] _ class.primitives[4] _ prim; <> x _ n + class.params[1]*N; y _ q + class.params[1]*Q; z _ l + class.params[1]*L; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; <> x _ n + class.params[2]*N; y _ q + class.params[2]*Q; z _ l + class.params[2]*L; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[2] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; x _ n + class.params[3]*N; y _ q + class.params[3]*Q; z _ l + class.params[3]*L; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[3] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; x _ n + class.params[4]*N; y _ q + class.params[4]*Q; z _ l + class.params[4]*L; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[4] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; }; ENDCASE => ERROR; }; SortHitRec: PROC [roots: Polynomial.ShortRealRootRec] = { <> temp: REAL; FOR i: NAT IN[0..roots.nRoots-1) DO FOR j: NAT IN[0..roots.nRoots-i-1) DO IF roots.realRoot[j] > roots.realRoot[j+1] THEN { temp _ roots.realRoot[j]; roots.realRoot[j] _ roots.realRoot[j+1]; roots.realRoot[j+1] _ temp;}; ENDLOOP; ENDLOOP; }; <> <> <> <> < roots[j+1] THEN {>> <> <> <> <> <> <<}; >> EvalTorus: PRIVATE PROC [RR, RRminusSS: REAL, degree: NAT, t: REAL, localRay: Ray] RETURNS [f: REAL] = { HH, x, y, z: REAL; p1,p2,p3,d1,d2,d3: REAL; basePt: Point3d; direction: Vector3d; [basePt, direction] _ SVRay.GetLocalRay[localRay]; p1 _ basePt[1]; p2 _ basePt[2]; p3 _ basePt[3]; d1 _ direction[1]; d2 _ direction[2]; d3 _ direction[3]; x _ p1 + t*d1; y _ p2 + t*d2; z _ p3 + t*d3; SELECT degree FROM 4=>{ xxPLUSzz: REAL; xxPLUSzz _ x*x + z*z; HH _ xxPLUSzz + y*y; f _ RRminusSS*RRminusSS + HH*(HH+2.0*RRminusSS) - 4.0*RR*xxPLUSzz; }; 3=>{ TwoHHdot: REAL; HH _ x*x + y*y + z*z; TwoHHdot _ 2.0*(x*d1 + y*d2 + z*d3); f _ TwoHHdot*(HH + 2.0*RRminusSS) + HH*TwoHHdot - 8.0*RR*(x*d1 + z*d3); }; 2=>{ HHdot, HHdotPrime: REAL; HH _ x*x + y*y + z*z; HHdot _ (x*d1 + y*d2 + z*d3); HHdotPrime _ (d1*d1+d2*d2+d3*d3); f _ 4.0*HHdotPrime*(HH + RRminusSS) + 8.0*(HHdot*HHdot - RR*(d1*d1+d3*d3)); }; 1=>{ f _ 24.0*(x*d1 + y*d2 + z*d3)*(d1*d1+d2*d2+d3*d3); }; ENDCASE => ERROR; }; TorusRoots: PRIVATE PROC [RR, RRminusSS: REAL, localRay: Ray, degree: NAT, bound: REAL] RETURNS [roots: Polynomial.ShortRealRootRec] = BEGIN roots.nRoots _ 0; IF degree=1 THEN BEGIN p: Point3d; d: Vector3d; [p,d] _ SVRay.GetLocalRay[localRay]; roots.nRoots _ 1; roots.realRoot[0] _ -(p[1]*d[1] + p[2]*d[2] + p[3]*d[3])/(d[1]*d[1] + d[2]*d[2] + d[3]*d[3]) END ELSE BEGIN BEGIN extrema: Polynomial.ShortRealRootRec _ TorusRoots[RR, RRminusSS, localRay, degree-1, bound]; x: REAL; <> <> <> <> <> IF extrema.nRoots>0 THEN BEGIN x _ TRootBetween[RR, RRminusSS, degree, -bound, extrema.realRoot[0], localRay]; IF x<=extrema.realRoot[0] THEN {roots.nRoots _ 1; roots.realRoot[0] _ x}; FOR i:NAT IN [0..extrema.nRoots-1) DO x _ TRootBetween[RR, RRminusSS, degree, extrema.realRoot[i], extrema.realRoot[i+1], localRay]; IF x<=extrema.realRoot[i+1] THEN {roots.realRoot[roots.nRoots] _ x; roots.nRoots _ roots.nRoots + 1}; ENDLOOP; x _ TRootBetween[RR, RRminusSS, degree, extrema.realRoot[extrema.nRoots-1],bound, localRay]; IF x<=bound THEN {roots.realRoot[roots.nRoots] _ x; roots.nRoots _ roots.nRoots + 1} END ELSE BEGIN x _ TRootBetween[RR, RRminusSS, degree, -bound, bound, localRay]; IF x<=bound THEN {roots.realRoot[0] _ x; roots.nRoots _ 1} END END END END; TRootBound: PROC [R, S: REAL, p: Point3d, d: Vector3d] RETURNS [REAL] = BEGIN -- finds a bound for the absolute values of the roots <> p1, dSup: REAL; p1 _ ABS[p[1]]+ABS[p[2]]+ABS[p[3]]; dSup _ MAX[ABS[d[1]], MAX[ABS[d[2]], ABS[d[3]]]]; RETURN[(p1 + R + S)/dSup]; END; TRootBetween: PROC [RR, RRminusSS: REAL, degree: NAT, x0, x1: REAL, localRay: Ray] RETURNS [x: REAL] = BEGIN y0: REAL _ EvalTorus[RR, RRminusSS, degree, x0, localRay]; y1: REAL _ EvalTorus[RR, RRminusSS, degree, x1, localRay]; <> <> xx: REAL; xx0: REAL _ IF y0yy0 THEN {xx0 _ x; yy0 _ y}}; IF y>0 THEN {IF y 0.0 AND y1 > 0.0) OR (y0 < 0.0 AND y1 < 0.0) THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval). < 0 THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval).>> 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 _ EvalTorus[RR, RRminusSS, degree, x0, localRay]; y1 _ EvalTorus[RR, RRminusSS, degree, x1, localRay]; THROUGH [0..500) DO y: REAL _ EvalTorus[RR, RRminusSS, degree, x _ (x0+x1)/2, localRay]; 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; PositiveTorusRoots: PRIVATE PROC [RR, RRminusSS: REAL, localRay: Ray, bound: REAL] RETURNS [roots: Polynomial.ShortRealRootRec] = { i: NAT _ 1; roots _ TorusRoots[RR, RRminusSS, localRay, 4, bound]; <> <> <> <> <> IF roots.nRoots = 0 THEN RETURN; WHILE i <= roots.nRoots DO IF roots.realRoot[i-1] <= 0 THEN { FOR j: NAT IN [i..roots.nRoots-1] DO roots.realRoot[j-1] _ roots.realRoot[j]; ENDLOOP; roots.nRoots _ roots.nRoots - 1; } ELSE {i _ i + 1}; ENDLOOP; }; InsideTorus: PRIVATE PROC [R, S: REAL, testPt: Point3d] RETURNS [BOOL] = { r, rMinusR: REAL; r _ RealFns.SqRt[testPt[1]*testPt[1] + testPt[3]*testPt[3]]; rMinusR _ r-R; RETURN[rMinusR*rMinusR + testPt[2]*testPt[2] <= S*S]; }; NewToroidCast: PUBLIC PROC [localRay: Ray, prim: Primitive, torusRec: REF ANY, surface: Surface] RETURNS [class: Classification] = { torusData: TorusRec; realRoots: Polynomial.ShortRealRootRec; -- optional B, R, S, RR, RRminusSS: REAL; -- big and little torus radii p: Point3d; d: Vector3d; bound: REAL; class _ SVCastRays.GetClassFromPool[]; torusData _ NARROW[torusRec]; R _ B _ torusData.bigRadius; S _ torusData.sectionRadius; RR _ R*R; RRminusSS _ RR - S*S; [p, d] _ SVRay.GetLocalRay[localRay]; bound _ TRootBound[R, S, p, d]; realRoots _ PositiveTorusRoots[RR, RRminusSS, localRay, bound]; SELECT realRoots.nRoots FROM <> 0 => {-- complete miss. class.count _ 0; class.classifs[1] _ FALSE; }; 1 => {-- a skim. Count as a miss. IF InsideTorus[R, S, p] THEN { x, y, z, r, OneMinusBoverR: REAL; class.count _ 1; class.classifs[1] _ TRUE; class.classifs[2] _ FALSE; class.params[1] _ realRoots.realRoot[0]; <> [x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[1]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; class.surfaces[1] _ surface; class.primitives[1] _ prim; } ELSE SVCastRays.MakeClassAMiss[class]; }; 2 => {-- cuts one cross section of the doughnut. Sort roots by size. x, y, z, r, OneMinusBoverR: REAL; class.count _ 2; IF realRoots.realRoot[0] < realRoots.realRoot[1] THEN {class.params[1] _ realRoots.realRoot[0]; class.params[2] _ realRoots.realRoot[1]} ELSE {class.params[2] _ realRoots.realRoot[0]; class.params[1] _ realRoots.realRoot[1]}; class.surfaces[1] _ surface; class.surfaces[2] _ surface; class.classifs[1] _ FALSE; class.classifs[2] _ TRUE; class.classifs[3] _ FALSE; class.primitives[1] _ class.primitives[2] _ prim; <> [x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[1]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; <> [x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[2]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[2] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; }; 3 => {-- cuts one cross section. Skims the other (or starts inside of the torus). For now, count as a miss. class.count _ 0; class.classifs[1] _ FALSE; }; 4 => {-- cuts both cross sections. Sort the roots. Pair them in order. x, y, z, r, OneMinusBoverR: REAL; <> SortHitRec[realRoots]; class.count _ 4; class.params[1] _ realRoots.realRoot[0]; class.params[2] _ realRoots.realRoot[1]; class.params[3] _ realRoots.realRoot[2]; class.params[4] _ realRoots.realRoot[3]; class.surfaces[1] _ class.surfaces[2] _ class.surfaces[3] _ class.surfaces[4] _ surface; class.classifs[1] _ FALSE; class.classifs[2] _ TRUE; class.classifs[3] _ FALSE; class.classifs[4] _ TRUE; class.classifs[5] _ FALSE; class.primitives[1] _ class.primitives[2] _ class.primitives[3] _ class.primitives[4] _ prim; <> [x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[1]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; <> [x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[2]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[2] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; [x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[3]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[3] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; [x, y, z] _ SVRay.EvaluateLocalRay2[localRay, class.params[4]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[4] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; }; ENDCASE => ERROR; }; 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.0 AND y1 > 0.0) OR (y0 < 0.0 AND y1 < 0.0) THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval). < 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*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; 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; 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; 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; 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; 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; Init: PROC = { globalPoly _ Polynomial.Quartic[[0,0,0,0,0]]; }; Init[]; END.