open-axiom repository from github
\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{src/algebra acplot.spad}
\author{Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\tableofcontents
\eject
\section{package REALSOLV RealSolvePackage}
<<package REALSOLV RealSolvePackage>>=
import PolynomialCategoryLifting
import Integer
import Float
import Symbol
import Fraction
import Polynomial
import FloatingRealPackage
)abbrev package REALSOLV RealSolvePackage
RealSolvePackage(): Exports == Implementation where
++ This package provides numerical solutions of systems of polynomial
++ equations for use in ACPLOT.
I ==> Integer
IE ==> IndexedExponents Symbol
L ==> List
NF ==> Float
P ==> Polynomial
RN ==> Fraction Integer
SE ==> Symbol
RFI ==> Fraction Polynomial Integer
LIFT ==> PolynomialCategoryLifting(IE,SE,RN,P RN,RFI)
SOLV ==> FloatingRealPackage Float
Exports ==> with
solve: (P RN,NF) -> L NF
++ solve(p,eps) finds the real zeroes of a
++ univariate rational polynomial p with precision eps.
solve: (P I,NF) -> L NF
++ solve(p,eps) finds the real zeroes of a univariate
++ integer polynomial p with precision eps.
realSolve: (L P I,L SE,NF) -> L L NF
++ realSolve(lp,lv,eps) = compute the list of the real
++ solutions of the list lp of polynomials with integer
++ coefficients with respect to the variables in lv,
++ with precision eps.
Implementation ==> add
prn2rfi: P RN -> RFI
prn2rfi p ==
map(#1 :: RFI,(numer(#1) :: RFI)/(denom(#1) :: RFI),p)$LIFT
pi2rfi: P I -> RFI
pi2rfi p == p :: RFI
solve(p:P RN,eps:NF) == realRoots(prn2rfi p,eps)$SOLV
solve(p:P I,eps:NF) ==
realRoots(p :: RFI,eps)$SOLV
realSolve(lp,lv,eps) ==
realRoots(map(pi2rfi,lp)$ListFunctions2(P I,RFI),lv,eps)$SOLV
@
\section{domain ACPLOT PlaneAlgebraicCurvePlot}
<<domain ACPLOT PlaneAlgebraicCurvePlot>>=
import PlottablePlaneCurveCategory
import Boolean
import Integer
import PositiveInteger
import DoubleFloat
import Float
import Symbol
import Polynomial
import UnivariatePolynomial
import SparseUnivariatePolynomial
import Point
import List
++ Plot a NON-SINGULAR plane algebraic curve p(x,y) = 0.
++ Author: Clifton J. Williamson
++ Date Created: Fall 1988
++ Date Last Updated: 27 April 1990
++ Keywords: algebraic curve, non-singular, plot
++ Examples:
++ References:
)abbrev domain ACPLOT PlaneAlgebraicCurvePlot
PlaneAlgebraicCurvePlot():Exports == Implementation where
B ==> Boolean
OUT ==> OutputForm
SE ==> Symbol
L ==> List
SEG ==> Segment
I ==> Integer
PI ==> PositiveInteger
RN ==> Fraction Integer
NF ==> Float
SF ==> DoubleFloat
P ==> Polynomial
UP ==> UnivariatePolynomial
SUP ==> SparseUnivariatePolynomial
FR ==> Factored
Pt ==> Point DoubleFloat
BoundaryPts ==> Record(left: L Pt,_
right: L Pt,_
bottom: L Pt,_
top: L Pt)
NewPtInfo ==> Record(newPt: Pt,_
type: String)
Corners ==> Record(minXVal: SF,_
maxXVal: SF,_
minYVal: SF,_
maxYVal: SF)
kinte ==> solve$RealSolvePackage()
rsolve ==> realSolve$RealSolvePackage()
Exports ==> PlottablePlaneCurveCategory with
makeSketch:(P I,SE,SE,SEG RN,SEG RN) -> %
++ makeSketch(p,x,y,a..b,c..d) creates an ACPLOT of the
++ curve \spad{p = 0} in the region {\em a <= x <= b, c <= y <= d}.
++ More specifically, 'makeSketch' plots a non-singular algebraic curve
++ \spad{p = 0} in an rectangular region {\em xMin <= x <= xMax},
++ {\em yMin <= y <= yMax}. The user inputs
++ \spad{makeSketch(p,x,y,xMin..xMax,yMin..yMax)}.
++ Here p is a polynomial in the variables x and y with
++ integer coefficients (p belongs to the domain
++ \spad{Polynomial Integer}). The case
++ where p is a polynomial in only one of the variables is
++ allowed. The variables x and y are input to specify the
++ the coordinate axes. The horizontal axis is the x-axis and
++ the vertical axis is the y-axis. The rational numbers
++ xMin,...,yMax specify the boundaries of the region in
++ which the curve is to be plotted.
refine:(%,SF) -> %
++ refine(p,x) \undocumented{}
Implementation ==> add
import PointPackage DoubleFloat
import Plot
import RealSolvePackage
singValBetween?:(SF,SF,L SF) -> B
segmentInfo:(SF -> SF,SF,SF,L SF,L SF,L SF,SF,SF) -> _
Record(seg:SEG SF,left: SF,lowerVals: L SF,upperVals:L SF)
swapCoords:Pt -> Pt
samePlottedPt?:(Pt,Pt) -> B
findPtOnList:(Pt,L Pt) -> Union(Pt,"failed")
makeCorners:(SF,SF,SF,SF) -> Corners
getXMin: Corners -> SF
getXMax: Corners -> SF
getYMin: Corners -> SF
getYMax: Corners -> SF
SFPolyToUPoly:P SF -> SUP SF
RNPolyToUPoly:P RN -> SUP RN
coerceCoefsToSFs:P I -> P SF
coerceCoefsToRNs:P I -> P RN
RNtoSF:RN -> SF
RNtoNF:RN -> NF
SFtoNF:SF -> NF
listPtsOnHorizBdry:(P RN,SE,RN,NF,NF) -> L Pt
listPtsOnVertBdry:(P RN,SE,RN,NF,NF) -> L Pt
listPtsInRect:(L L NF,NF,NF,NF,NF) -> L Pt
ptsSuchThat?:(L L NF,L NF -> B) -> B
inRect?:(L NF,NF,NF,NF,NF) -> B
onHorzSeg?:(L NF,NF,NF,NF) -> B
onVertSeg?:(L NF,NF,NF,NF) -> B
newX:(L L NF,L L NF,NF,NF,NF,RN,RN) -> RN
newY:(L L NF,L L NF,NF,NF,NF,RN,RN) -> RN
makeOneVarSketch:(P I,SE,SE,RN,RN,RN,RN,SE) -> %
makeLineSketch:(P I,SE,SE,RN,RN,RN,RN) -> %
makeRatFcnSketch:(P I,SE,SE,RN,RN,RN,RN,SE) -> %
makeGeneralSketch:(P I,SE,SE,RN,RN,RN,RN) -> %
traceBranches:(P SF,P SF,P SF,SE,SE,Corners,SF,SF,PI,_
L Pt,BoundaryPts) -> L L Pt
dummyFirstPt:(Pt,P SF,P SF,SE,SE,L Pt,L Pt,L Pt,L Pt) -> Pt
listPtsOnSegment:(P SF,P SF,P SF,SE,SE,Pt,Pt,Corners,_
SF,SF,PI,L Pt,L Pt) -> L L Pt
listPtsOnLoop:(P SF,P SF,P SF,SE,SE,Pt,Corners,_
SF,SF,PI,L Pt,L Pt) -> L L Pt
computeNextPt:(P SF,P SF,P SF,SE,SE,Pt,Pt,Corners,_
SF,SF,PI,L Pt,L Pt) -> NewPtInfo
newtonApprox:(SUP SF, SF, SF, PI) -> Union(SF, "failed")
--% representation
Rep := Record(poly : P I,_
xVar : SE,_
yVar : SE,_
minXVal : RN,_
maxXVal : RN,_
minYVal : RN,_
maxYVal : RN,_
bdryPts : BoundaryPts,_
hTanPts : L Pt,_
vTanPts : L Pt,_
branches: L L Pt)
--% global constants
EPSILON : NF := 0.000001 -- precision to which realSolve finds roots
PLOTERR : SF := float(1,-3,10)
-- maximum allowable difference in each coordinate when
-- determining if 2 plotted points are equal
--% global flags
NADA : String := "nothing in particular"
BDRY : String := "boundary point"
CRIT : String := "critical point"
BOTTOM : String := "bottom"
TOP : String := "top"
--% hacks
NFtoSF: NF -> SF
NFtoSF x == 0 + convert(x)$NF
--% points
makePt: (SF,SF) -> Pt
makePt(xx,yy) == point(l : L SF := [xx,yy])
swapCoords(pt) == makePt(yCoord pt,xCoord pt)
samePlottedPt?(p0,p1) ==
-- determines if p1 lies in a square with side 2 PLOTERR
-- centered at p0
x0 := xCoord p0; y0 := yCoord p0
x1 := xCoord p1; y1 := yCoord p1
(abs(x1-x0) < PLOTERR) and (abs(y1-y0) < PLOTERR)
findPtOnList(pt,pointList) ==
for point in pointList repeat
samePlottedPt?(pt,point) => return point
"failed"
--% corners
makeCorners(xMinSF,xMaxSF,yMinSF,yMaxSF) ==
[xMinSF,xMaxSF,yMinSF,yMaxSF]
getXMin(corners) == corners.minXVal
getXMax(corners) == corners.maxXVal
getYMin(corners) == corners.minYVal
getYMax(corners) == corners.maxYVal
--% coercions
SFPolyToUPoly(p) ==
-- 'p' is of type Polynomial, but has only one variable
zero? p => 0
monomial(leadingCoefficient p,totalDegree p) +
SFPolyToUPoly(reductum p)
RNPolyToUPoly(p) ==
-- 'p' is of type Polynomial, but has only one variable
zero? p => 0
monomial(leadingCoefficient p,totalDegree p) +
RNPolyToUPoly(reductum p)
coerceCoefsToSFs(p) ==
-- coefficients of 'p' are coerced to be DoubleFloat's
map(coerce,p)$PolynomialFunctions2(I,SF)
coerceCoefsToRNs(p) ==
-- coefficients of 'p' are coerced to be DoubleFloat's
map(coerce,p)$PolynomialFunctions2(I,RN)
RNtoSF(r) == coerce(r)@SF
RNtoNF(r) == coerce(r)@NF
SFtoNF(x) == convert(x)@NF
--% computation of special points
listPtsOnHorizBdry(pRN,y,y0,xMinNF,xMaxNF) ==
-- strict inequality here: corners on vertical boundary
pointList : L Pt := nil()
ySF := RNtoSF(y0)
f := eval(pRN,y,y0)
roots : L NF := kinte(f,EPSILON)
for root in roots repeat
if (xMinNF < root) and (root < xMaxNF) then
pointList := cons(makePt(NFtoSF root, ySF), pointList)
pointList
listPtsOnVertBdry(pRN,x,x0,yMinNF,yMaxNF) ==
pointList : L Pt := nil()
xSF := RNtoSF(x0)
f := eval(pRN,x,x0)
roots : L NF := kinte(f,EPSILON)
for root in roots repeat
if (yMinNF <= root) and (root <= yMaxNF) then
pointList := cons(makePt(xSF, NFtoSF root), pointList)
pointList
listPtsInRect(points,xMin,xMax,yMin,yMax) ==
pointList : L Pt := nil()
for point in points repeat
xx := first point; yy := second point
if (xMin<=xx) and (xx<=xMax) and (yMin<=yy) and (yy<=yMax) then
pointList := cons(makePt(NFtoSF xx,NFtoSF yy),pointList)
pointList
ptsSuchThat?(points,pred) ==
for point in points repeat
if pred point then return true
false
inRect?(point,xMinNF,xMaxNF,yMinNF,yMaxNF) ==
xx := first point; yy := second point
xMinNF <= xx and xx <= xMaxNF and yMinNF <= yy and yy <= yMaxNF
onHorzSeg?(point,xMinNF,xMaxNF,yNF) ==
xx := first point; yy := second point
yy = yNF and xMinNF <= xx and xx <= xMaxNF
onVertSeg?(point,yMinNF,yMaxNF,xNF) ==
xx := first point; yy := second point
xx = xNF and yMinNF <= yy and yy <= yMaxNF
newX(vtanPts,singPts,yMinNF,yMaxNF,xNF,xRN,horizInc) ==
xNewNF := xNF + RNtoNF horizInc
xRtNF := max(xNF,xNewNF); xLftNF := min(xNF,xNewNF)
-- ptsSuchThat?(singPts,inRect?(#1,xLftNF,xRtNF,yMinNF,yMaxNF)) =>
foo : L NF -> B := inRect?(#1,xLftNF,xRtNF,yMinNF,yMaxNF)
ptsSuchThat?(singPts,foo) =>
newX(vtanPts,singPts,yMinNF,yMaxNF,xNF,xRN,horizInc/2::RN)
-- ptsSuchThat?(vtanPts,onVertSeg?(#1,yMinNF,yMaxNF,xNewNF)) =>
goo : L NF -> B := onVertSeg?(#1,yMinNF,yMaxNF,xNewNF)
ptsSuchThat?(vtanPts,goo) =>
newX(vtanPts,singPts,yMinNF,yMaxNF,xNF,xRN,horizInc/2::RN)
xRN + horizInc
newY(htanPts,singPts,xMinNF,xMaxNF,yNF,yRN,vertInc) ==
yNewNF := yNF + RNtoNF vertInc
yTopNF := max(yNF,yNewNF); yBotNF := min(yNF,yNewNF)
-- ptsSuchThat?(singPts,inRect?(#1,xMinNF,xMaxNF,yBotNF,yTopNF)) =>
foo : L NF -> B := inRect?(#1,xMinNF,xMaxNF,yBotNF,yTopNF)
ptsSuchThat?(singPts,foo) =>
newY(htanPts,singPts,xMinNF,xMaxNF,yNF,yRN,vertInc/2::RN)
-- ptsSuchThat?(htanPts,onHorzSeg?(#1,xMinNF,xMaxNF,yNewNF)) =>
goo : L NF -> B := onHorzSeg?(#1,xMinNF,xMaxNF,yNewNF)
ptsSuchThat?(htanPts,goo) =>
newY(htanPts,singPts,xMinNF,xMaxNF,yNF,yRN,vertInc/2::RN)
yRN + vertInc
--% creation of sketches
makeSketch(p,x,y,xRange,yRange) ==
xMin := lo xRange; xMax := hi xRange
yMin := lo yRange; yMax := hi yRange
-- test input for consistency
xMax <= xMin =>
error "makeSketch: bad range for first variable"
yMax <= yMin =>
error "makeSketch: bad range for second variable"
varList := variables p
# varList > 2 =>
error "makeSketch: polynomial in more than 2 variables"
# varList = 0 =>
error "makeSketch: constant polynomial"
-- polynomial in 1 variable
# varList = 1 =>
(not member?(x,varList)) and (not member?(y,varList)) =>
error "makeSketch: bad variables"
makeOneVarSketch(p,x,y,xMin,xMax,yMin,yMax,first varList)
-- polynomial in 2 variables
(not member?(x,varList)) or (not member?(y,varList)) =>
error "makeSketch: bad variables"
totalDegree p = 1 =>
makeLineSketch(p,x,y,xMin,xMax,yMin,yMax)
-- polynomial is linear in one variable
-- y is a rational function of x
degree(p,y) = 1 =>
makeRatFcnSketch(p,x,y,xMin,xMax,yMin,yMax,y)
-- x is a rational function of y
degree(p,x) = 1 =>
makeRatFcnSketch(p,x,y,xMin,xMax,yMin,yMax,x)
-- the general case
makeGeneralSketch(p,x,y,xMin,xMax,yMin,yMax)
--% special cases
makeOneVarSketch(p,x,y,xMin,xMax,yMin,yMax,var) ==
-- the case where 'p' is a polynomial in only one variable
-- the graph consists of horizontal or vertical lines
if var = x then
minVal := RNtoNF xMin
maxVal := RNtoNF xMax
else
minVal := RNtoNF yMin
maxVal := RNtoNF yMax
lf : L Pt := nil(); rt : L Pt := nil()
bt : L Pt := nil(); tp : L Pt := nil()
htans : L Pt := nil(); vtans : L Pt := nil()
bran : L L Pt := nil()
roots := kinte(p,EPSILON)
sketchRoots : L SF := nil()
for root in roots repeat
if (minVal <= root) and (root <= maxVal) then
sketchRoots := cons(NFtoSF root,sketchRoots)
null sketchRoots =>
[p,x,y,xMin,xMax,yMin,yMax,[lf,rt,bt,tp],htans,vtans,bran]
if var = x then
yMinSF := RNtoSF yMin; yMaxSF := RNtoSF yMax
for rootSF in sketchRoots repeat
tp := cons(pt1 := makePt(rootSF,yMaxSF),tp)
bt := cons(pt2 := makePt(rootSF,yMinSF),bt)
branch : L Pt := [pt1,pt2]
bran := cons(branch,bran)
else
xMinSF := RNtoSF xMin; xMaxSF := RNtoSF xMax
for rootSF in sketchRoots repeat
rt := cons(pt1 := makePt(xMaxSF,rootSF),rt)
lf := cons(pt2 := makePt(xMinSF,rootSF),lf)
branch : L Pt := [pt1,pt2]
bran := cons(branch,bran)
[p,x,y,xMin,xMax,yMin,yMax,[lf,rt,bt,tp],htans,vtans,bran]
makeLineSketch(p,x,y,xMin,xMax,yMin,yMax) ==
-- the case where p(x,y) = a x + b y + c with a ~= 0, b ~= 0
-- this is a line which is neither vertical nor horizontal
xMinSF := RNtoSF xMin; xMaxSF := RNtoSF xMax
yMinSF := RNtoSF yMin; yMaxSF := RNtoSF yMax
-- determine the coefficients a, b, and c
a := ground(coefficient(p,x,1)) :: SF
b := ground(coefficient(p,y,1)) :: SF
c := ground(coefficient(coefficient(p,x,0),y,0)) :: SF
lf : L Pt := nil(); rt : L Pt := nil()
bt : L Pt := nil(); tp : L Pt := nil()
htans : L Pt := nil(); vtans : L Pt := nil()
branch : L Pt := nil(); bran : L L Pt := nil()
-- compute x coordinate of point on line with y = yMin
xBottom := (- b*yMinSF - c)/a
-- compute x coordinate of point on line with y = yMax
xTop := (- b*yMaxSF - c)/a
-- compute y coordinate of point on line with x = xMin
yLeft := (- a*xMinSF - c)/b
-- compute y coordinate of point on line with x = xMax
yRight := (- a*xMaxSF - c)/b
-- determine which of the above 4 points are in the region
-- to be plotted and list them as a branch
if (xMinSF < xBottom) and (xBottom < xMaxSF) then
bt := cons(pt := makePt(xBottom,yMinSF),bt)
branch := cons(pt,branch)
if (xMinSF < xTop) and (xTop < xMaxSF) then
tp := cons(pt := makePt(xTop,yMaxSF),tp)
branch := cons(pt,branch)
if (yMinSF <= yLeft) and (yLeft <= yMaxSF) then
lf := cons(pt := makePt(xMinSF,yLeft),lf)
branch := cons(pt,branch)
if (yMinSF <= yRight) and (yRight <= yMaxSF) then
rt := cons(pt := makePt(xMaxSF,yRight),rt)
branch := cons(pt,branch)
bran := cons(branch,bran)
[p,x,y,xMin,xMax,yMin,yMax,[lf,rt,bt,tp],htans,vtans,bran]
singValBetween?(xCurrent,xNext,xSingList) ==
for xVal in xSingList repeat
(xCurrent < xVal) and (xVal < xNext) => return true
false
segmentInfo(f,lo,hi,botList,topList,singList,minSF,maxSF) ==
repeat
-- 'current' is the smallest element of 'topList' and 'botList'
-- 'currentFrom' records the list from which it was taken
if null topList then
if null botList then
return [segment(lo,hi),hi,nil(),nil()]
else
current := first botList
botList := rest botList
currentFrom := BOTTOM
else
if null botList then
current := first topList
topList := rest topList
currentFrom := TOP
else
bot := first botList
top := first topList
if bot < top then
current := bot
botList := rest botList
currentFrom := BOTTOM
else
current := top
topList := rest topList
currentFrom := TOP
-- 'nxt' is the next smallest element of 'topList'
-- and 'botList'
-- 'nextFrom' records the list from which it was taken
if null topList then
if null botList then
return [segment(lo,hi),hi,nil(),nil()]
else
nxt := first botList
botList := rest botList
nextFrom := BOTTOM
else
if null botList then
nxt := first topList
topList := rest topList
nextFrom := TOP
else
bot := first botList
top := first topList
if bot < top then
nxt := bot
botList := rest botList
nextFrom := BOTTOM
else
nxt := top
topList := rest topList
nextFrom := TOP
if currentFrom = nextFrom then
if singValBetween?(current,nxt,singList) then
return [segment(lo,current),nxt,botList,topList]
else
val := f((nxt - current)/2::SF)
if (val <= minSF) or (val >= maxSF) then
return [segment(lo,current),nxt,botList,topList]
else
if singValBetween?(current,nxt,singList) then
return [segment(lo,current),nxt,botList,topList]
makeRatFcnSketch(p,x,y,xMin,xMax,yMin,yMax,depVar) ==
-- the case where p(x,y) is linear in x or y
-- Thus, one variable is a rational function of the other.
-- Therefore, we may use the 2-dimensional function plotting
-- package. The only problem is determining the intervals on
-- on which the function is to be plotted.
--!! corners: e.g. upper left corner is on graph with y' > 0
factoredP := p :: FR P I
numberOfFactors(factoredP) > 1 =>
error "reducible polynomial" --!! sketch each factor
dpdx := differentiate(p,x)
dpdy := differentiate(p,y)
pRN := coerceCoefsToRNs p
xMinSF := RNtoSF xMin; xMaxSF := RNtoSF xMax
yMinSF := RNtoSF yMin; yMaxSF := RNtoSF yMax
xMinNF := RNtoNF xMin; xMaxNF := RNtoNF xMax
yMinNF := RNtoNF yMin; yMaxNF := RNtoNF yMax
-- 'p' is of degree 1 in the variable 'depVar'.
-- Thus, 'depVar' is a rational function of the other variable.
num := -coefficient(p,depVar,0)
den := coefficient(p,depVar,1)
numUPolySF := SFPolyToUPoly(coerceCoefsToSFs(num))
denUPolySF := SFPolyToUPoly(coerceCoefsToSFs(den))
-- this is the rational function
f : SF -> SF := elt(numUPolySF,#1)/elt(denUPolySF,#1)
-- values of the dependent and independent variables
if depVar = x then
indVarMin := yMin; indVarMax := yMax
indVarMinNF := yMinNF; indVarMaxNF := yMaxNF
indVarMinSF := yMinSF; indVarMaxSF := yMaxSF
depVarMin := xMin; depVarMax := xMax
depVarMinSF := xMinSF; depVarMaxSF := xMaxSF
else
indVarMin := xMin; indVarMax := xMax
indVarMinNF := xMinNF; indVarMaxNF := xMaxNF
indVarMinSF := xMinSF; indVarMaxSF := xMaxSF
depVarMin := yMin; depVarMax := yMax
depVarMinSF := yMinSF; depVarMaxSF := yMaxSF
-- Create lists of critical points.
htanPts := rsolve([p,dpdx],[x,y],EPSILON)
vtanPts := rsolve([p,dpdy],[x,y],EPSILON)
htans := listPtsInRect(htanPts,xMinNF,xMaxNF,yMinNF,yMaxNF)
vtans := listPtsInRect(vtanPts,xMinNF,xMaxNF,yMinNF,yMaxNF)
-- Create lists which will contain boundary points.
lf : L Pt := nil(); rt : L Pt := nil()
bt : L Pt := nil(); tp : L Pt := nil()
-- Determine values of the independent variable at the which
-- the rational function has a pole as well as the values of
-- the independent variable for which there is a point on the
-- upper or lower boundary.
singList : L SF :=
roots : L NF := kinte(den,EPSILON)
outList : L SF := nil()
for root in roots repeat
if (indVarMinNF < root) and (root < indVarMaxNF) then
outList := cons(NFtoSF root,outList)
sort(#1 < #2,outList)
topList : L SF :=
roots : L NF := kinte(eval(pRN,depVar,depVarMax),EPSILON)
outList : L SF := nil()
for root in roots repeat
if (indVarMinNF < root) and (root < indVarMaxNF) then
outList := cons(NFtoSF root,outList)
sort(#1 < #2,outList)
botList : L SF :=
roots : L NF := kinte(eval(pRN,depVar,depVarMin),EPSILON)
outList : L SF := nil()
for root in roots repeat
if (indVarMinNF < root) and (root < indVarMaxNF) then
outList := cons(NFtoSF root,outList)
sort(#1 < #2,outList)
-- We wish to determine if the graph has points on the 'left'
-- and 'right' boundaries, so we compute the value of the
-- rational function at the lefthand and righthand values of
-- the dependent variable. If the function has a singularity
-- on the left or right boundary, then 'leftVal' or 'rightVal'
-- is given a dummy valuewhich will convince the program that
-- there is no point on the left or right boundary.
denUPolyRN := RNPolyToUPoly(coerceCoefsToRNs(den))
if elt(denUPolyRN,indVarMin) = 0$RN then
leftVal := depVarMinSF - (abs(depVarMinSF) + 1$SF)
else
leftVal := f(indVarMinSF)
if elt(denUPolyRN,indVarMax) = 0$RN then
rightVal := depVarMinSF - (abs(depVarMinSF) + 1$SF)
else
rightVal := f(indVarMaxSF)
-- Now put boundary points on the appropriate lists.
if depVar = x then
if (xMinSF < leftVal) and (leftVal < xMaxSF) then
bt := cons(makePt(leftVal,yMinSF),bt)
if (xMinSF < rightVal) and (rightVal < xMaxSF) then
tp := cons(makePt(rightVal,yMaxSF),tp)
for val in botList repeat
lf := cons(makePt(xMinSF,val),lf)
for val in topList repeat
rt := cons(makePt(xMaxSF,val),rt)
else
if (yMinSF < leftVal) and (leftVal < yMaxSF) then
lf := cons(makePt(xMinSF,leftVal),lf)
if (yMinSF < rightVal) and (rightVal < yMaxSF) then
rt := cons(makePt(xMaxSF,rightVal),rt)
for val in botList repeat
bt := cons(makePt(val,yMinSF),bt)
for val in topList repeat
tp := cons(makePt(val,yMaxSF),tp)
bran : L L Pt := nil()
-- Determine segments on which the rational function is to
-- be plotted.
if (depVarMinSF < leftVal) and (leftVal < depVarMaxSF) then
lo := indVarMinSF
else
if null topList then
if null botList then
return [p,x,y,xMin,xMax,yMin,yMax,[lf,rt,bt,tp],_
htans,vtans,bran]
else
lo := first botList
botList := rest botList
else
if null botList then
lo := first topList
topList := rest topList
else
bot := first botList
top := first topList
if bot < top then
lo := bot
botList := rest botList
else
lo := top
topList := rest topList
hi := 0$SF -- @#$%^&* compiler
if (depVarMinSF < rightVal) and (rightVal < depVarMaxSF) then
hi := indVarMaxSF
else
if null topList then
if null botList then
error "makeRatFcnSketch: plot domain"
else
hi := last botList
botList := remove(hi,botList)
else
if null botList then
hi := last topList
topList := remove(hi,topList)
else
bot := last botList
top := last topList
if bot > top then
hi := bot
botList := remove(hi,botList)
else
hi := top
topList := remove(hi,topList)
if (depVar = x) then
(minSF := xMinSF; maxSF := xMaxSF)
else
(minSF := yMinSF; maxSF := yMaxSF)
segList : L SEG SF := nil()
repeat
segInfo := segmentInfo(f,lo,hi,botList,topList,singList,_
minSF,maxSF)
segList := cons(segInfo.seg,segList)
lo := segInfo.left
botList := segInfo.lowerVals
topList := segInfo.upperVals
if lo = hi then break
for segment in segList repeat
RFPlot : Plot := plot(f,segment)
curve := first(listBranches(RFPlot))
if depVar = y then
bran := cons(curve,bran)
else
bran := cons(map(swapCoords,curve),bran)
[p,x,y,xMin,xMax,yMin,yMax,[lf,rt,bt,tp],htans,vtans,bran]
--% the general case
makeGeneralSketch(pol,x,y,xMin,xMax,yMin,yMax) ==
--!! corners of region should not be on curve
--!! enlarge region if necessary
factoredPol := pol :: FR P I
numberOfFactors(factoredPol) > 1 =>
error "reducible polynomial" --!! sketch each factor
p := nthFactor(factoredPol,1)
dpdx := differentiate(p,x); dpdy := differentiate(p,y)
xMinNF := RNtoNF xMin; xMaxNF := RNtoNF xMax
yMinNF := RNtoNF yMin; yMaxNF := RNtoNF yMax
-- compute singular points; error if singularities in region
singPts := rsolve([p,dpdx,dpdy],[x,y],EPSILON)
-- ptsSuchThat?(singPts,inRect?(#1,xMinNF,xMaxNF,yMinNF,yMaxNF)) =>
foo : L NF -> B := inRect?(#1,xMinNF,xMaxNF,yMinNF,yMaxNF)
ptsSuchThat?(singPts,foo) =>
error "singular pts in region of sketch"
-- compute critical points
htanPts := rsolve([p,dpdx],[x,y],EPSILON)
vtanPts := rsolve([p,dpdy],[x,y],EPSILON)
critPts := append(htanPts,vtanPts)
-- if there are critical points on the boundary, then enlarge
-- the region, but be sure that the new region does not contain
-- any singular points
hInc : RN := (1/20) * (xMax - xMin)
vInc : RN := (1/20) * (yMax - yMin)
-- if ptsSuchThat?(critPts,onVertSeg?(#1,yMinNF,yMaxNF,xMinNF)) then
foo : L NF -> B := onVertSeg?(#1,yMinNF,yMaxNF,xMinNF)
if ptsSuchThat?(critPts,foo) then
xMin := newX(critPts,singPts,yMinNF,yMaxNF,xMinNF,xMin,-hInc)
xMinNF := RNtoNF xMin
-- if ptsSuchThat?(critPts,onVertSeg?(#1,yMinNF,yMaxNF,xMaxNF)) then
foo : L NF -> B := onVertSeg?(#1,yMinNF,yMaxNF,xMaxNF)
if ptsSuchThat?(critPts,foo) then
xMax := newX(critPts,singPts,yMinNF,yMaxNF,xMaxNF,xMax,hInc)
xMaxNF := RNtoNF xMax
-- if ptsSuchThat?(critPts,onHorzSeg?(#1,xMinNF,xMaxNF,yMinNF)) then
foo : L NF -> B := onHorzSeg?(#1,xMinNF,xMaxNF,yMinNF)
if ptsSuchThat?(critPts,foo) then
yMin := newY(critPts,singPts,xMinNF,xMaxNF,yMinNF,yMin,-vInc)
yMinNF := RNtoNF yMin
-- if ptsSuchThat?(critPts,onHorzSeg?(#1,xMinNF,xMaxNF,yMaxNF)) then
foo : L NF -> B := onHorzSeg?(#1,xMinNF,xMaxNF,yMaxNF)
if ptsSuchThat?(critPts,foo) then
yMax := newY(critPts,singPts,xMinNF,xMaxNF,yMaxNF,yMax,vInc)
yMaxNF := RNtoNF yMax
htans := listPtsInRect(htanPts,xMinNF,xMaxNF,yMinNF,yMaxNF)
vtans := listPtsInRect(vtanPts,xMinNF,xMaxNF,yMinNF,yMaxNF)
crits := append(htans,vtans)
-- conversions to DoubleFloats
xMinSF := RNtoSF xMin; xMaxSF := RNtoSF xMax
yMinSF := RNtoSF yMin; yMaxSF := RNtoSF yMax
corners := makeCorners(xMinSF,xMaxSF,yMinSF,yMaxSF)
pSF := coerceCoefsToSFs p
dpdxSF := coerceCoefsToSFs dpdx
dpdySF := coerceCoefsToSFs dpdy
delta := min((xMaxSF - xMinSF)/25,(yMaxSF - yMinSF)/25)
err := min(delta/100,PLOTERR/100)
bound : PI := 10
-- compute points on the boundary
pRN := coerceCoefsToRNs(p)
lf : L Pt := listPtsOnVertBdry(pRN,x,xMin,yMinNF,yMaxNF)
rt : L Pt := listPtsOnVertBdry(pRN,x,xMax,yMinNF,yMaxNF)
bt : L Pt := listPtsOnHorizBdry(pRN,y,yMin,xMinNF,xMaxNF)
tp : L Pt := listPtsOnHorizBdry(pRN,y,yMax,xMinNF,xMaxNF)
bdPts : BoundaryPts := [lf,rt,bt,tp]
bran := traceBranches(pSF,dpdxSF,dpdySF,x,y,corners,delta,err,_
bound,crits,bdPts)
[p,x,y,xMin,xMax,yMin,yMax,bdPts,htans,vtans,bran]
refine(plot,stepFraction) ==
p := plot.poly; x := plot.xVar; y := plot.yVar
dpdx := differentiate(p,x); dpdy := differentiate(p,y)
pSF := coerceCoefsToSFs p
dpdxSF := coerceCoefsToSFs dpdx
dpdySF := coerceCoefsToSFs dpdy
xMin := plot.minXVal; xMax := plot.maxXVal
yMin := plot.minYVal; yMax := plot.maxYVal
xMinSF := RNtoSF xMin; xMaxSF := RNtoSF xMax
yMinSF := RNtoSF yMin; yMaxSF := RNtoSF yMax
corners := makeCorners(xMinSF,xMaxSF,yMinSF,yMaxSF)
pSF := coerceCoefsToSFs p
dpdxSF := coerceCoefsToSFs dpdx
dpdySF := coerceCoefsToSFs dpdy
delta :=
stepFraction * min((xMaxSF - xMinSF)/25,(yMaxSF - yMinSF)/25)
err := min(delta/100,PLOTERR/100)
bound : PI := 10
crits := append(plot.hTanPts,plot.vTanPts)
bdPts := plot.bdryPts
bran := traceBranches(pSF,dpdxSF,dpdySF,x,y,corners,delta,err,_
bound,crits,bdPts)
htans := plot.hTanPts; vtans := plot.vTanPts
[p,x,y,xMin,xMax,yMin,yMax,bdPts,htans,vtans,bran]
traceBranches(pSF,dpdxSF,dpdySF,x,y,corners,delta,err,bound,_
crits,bdPts) ==
-- for boundary points, trace curve from boundary to boundary
-- add the branch to the list of branches
-- update list of boundary points by deleting first and last
-- points on this branch
-- update list of critical points by deleting any critical
-- points which were plotted
lf := bdPts.left; rt := bdPts.right
tp := bdPts.top ; bt := bdPts.bottom
bdry := append(append(lf,rt),append(bt,tp))
bran : L L Pt := nil()
while not null bdry repeat
pt := first bdry
p0 := dummyFirstPt(pt,dpdxSF,dpdySF,x,y,lf,rt,bt,tp)
segInfo := listPtsOnSegment(pSF,dpdxSF,dpdySF,x,y,p0,pt,_
corners,delta,err,bound,crits,bdry)
bran := cons(first segInfo,bran)
crits := second segInfo
bdry := third segInfo
-- trace loops beginning and ending with critical points
-- add the branch to the list of branches
-- update list of critical points by deleting any critical
-- points which were plotted
while not null crits repeat
pt := first crits
segInfo := listPtsOnLoop(pSF,dpdxSF,dpdySF,x,y,pt,_
corners,delta,err,bound,crits,bdry)
bran := cons(first segInfo,bran)
crits := second segInfo
bran
dummyFirstPt(p1,dpdxSF,dpdySF,x,y,lf,rt,bt,tp) ==
-- The function 'computeNextPt' requires 2 points, p0 and p1.
-- When computing the second point on a branch which starts
-- on the boundary, we use the boundary point as p1 and the
-- 'dummy' point returned by this function as p0.
x1 := xCoord p1; y1 := yCoord p1
zero := 0$SF; one := 1$SF
px := ground(eval(dpdxSF,[x,y],[x1,y1]))
py := ground(eval(dpdySF,[x,y],[x1,y1]))
if px * py < zero then -- positive slope at p1
member?(p1,lf) or member?(p1,bt) =>
makePt(x1 - one,y1 - one)
makePt(x1 + one,y1 + one)
else
member?(p1,lf) or member?(p1,tp) =>
makePt(x1 - one,y1 + one)
makePt(x1 + one,y1 - one)
listPtsOnSegment(pSF,dpdxSF,dpdySF,x,y,p0,p1,corners,_
delta,err,bound,crits,bdry) ==
-- p1 is a boundary point; p0 is a 'dummy' point
bdry := remove(p1,bdry)
pointList : L Pt := [p1]
ptInfo := computeNextPt(pSF,dpdxSF,dpdySF,x,y,p0,p1,corners,_
delta,err,bound,crits,bdry)
p2 := ptInfo.newPt
ptInfo.type = BDRY =>
bdry := remove(p2,bdry)
pointList := cons(p2,pointList)
[pointList,crits,bdry]
if ptInfo.type = CRIT then crits := remove(p2,crits)
pointList := cons(p2,pointList)
repeat
pt0 := second pointList; pt1 := first pointList
ptInfo := computeNextPt(pSF,dpdxSF,dpdySF,x,y,pt0,pt1,corners,_
delta,err,bound,crits,bdry)
p2 := ptInfo.newPt
ptInfo.type = BDRY =>
bdry := remove(p2,bdry)
pointList := cons(p2,pointList)
return [pointList,crits,bdry]
if ptInfo.type = CRIT then crits := remove(p2,crits)
pointList := cons(p2,pointList)
--!! delete next line (compiler bug)
[pointList,crits,bdry]
listPtsOnLoop(pSF,dpdxSF,dpdySF,x,y,p1,corners,_
delta,err,bound,crits,bdry) ==
x1 := xCoord p1; y1 := yCoord p1
px := ground(eval(dpdxSF,[x,y],[x1,y1]))
py := ground(eval(dpdySF,[x,y],[x1,y1]))
p0 := makePt(x1 - 1$SF,y1 - 1$SF)
pointList : L Pt := [p1]
ptInfo := computeNextPt(pSF,dpdxSF,dpdySF,x,y,p0,p1,corners,_
delta,err,bound,crits,bdry)
p2 := ptInfo.newPt
ptInfo.type = BDRY =>
error "boundary reached while on loop"
if ptInfo.type = CRIT then
p1 = p2 =>
error "first and second points on loop are identical"
crits := remove(p2,crits)
pointList := cons(p2,pointList)
repeat
pt0 := second pointList; pt1 := first pointList
ptInfo := computeNextPt(pSF,dpdxSF,dpdySF,x,y,pt0,pt1,corners,_
delta,err,bound,crits,bdry)
p2 := ptInfo.newPt
ptInfo.type = BDRY =>
error "boundary reached while on loop"
if ptInfo.type = CRIT then
crits := remove(p2,crits)
p1 = p2 =>
pointList := cons(p2,pointList)
return [pointList,crits,bdry]
pointList := cons(p2,pointList)
--!! delete next line (compiler bug)
[pointList,crits,bdry]
computeNextPt(pSF,dpdxSF,dpdySF,x,y,p0,p1,corners,_
delta,err,bound,crits,bdry) ==
-- p0=(x0,y0) and p1=(x1,y1) are the last two points on the curve.
-- The function computes the next point on the curve.
-- The function determines if the next point is a critical point
-- or a boundary point.
-- The function returns a record of the form
-- Record(newPt:Pt,type:String).
-- If the new point is a boundary point, then 'type' is
-- "boundary point" and 'newPt' is a boundary point to be
-- deleted from the list of boundary points yet to be plotted.
-- Similarly, if the new point is a critical point, then 'type' is
-- "critical point" and 'newPt' is a critical point to be
-- deleted from the list of critical points yet to be plotted.
-- If the new point is neither a critical point nor a boundary
-- point, then 'type' is "nothing in particular".
xMinSF := getXMin corners; xMaxSF := getXMax corners
yMinSF := getYMin corners; yMaxSF := getYMax corners
x0 := xCoord p0; y0 := yCoord p0
x1 := xCoord p1; y1 := yCoord p1
px := ground(eval(dpdxSF,[x,y],[x1,y1]))
py := ground(eval(dpdySF,[x,y],[x1,y1]))
-- let m be the slope of the tangent line at p1
-- if |m| < 1, we will increment the x-coordinate by delta
-- (indicated by 'incVar = x'), find an approximate
-- y-coordinate using the tangent line, then find the actual
-- y-coordinate using a Newton iteration
if abs(py) > abs(px) then
incVar0 := incVar := x
deltaX := (if x1 > x0 then delta else -delta)
x2Approx := x1 + deltaX
y2Approx := y1 + (-px/py)*deltaX
-- if |m| >= 1, we interchange the roles of the x- and y-
-- coordinates
else
incVar0 := incVar := y
deltaY := (if y1 > y0 then delta else -delta)
x2Approx := x1 + (-py/px)*deltaY
y2Approx := y1 + deltaY
lookingFor := NADA
-- See if (x2Approx,y2Approx) is out of bounds.
-- If so, find where the line segment connecting (x1,y1) and
-- (x2Approx,y2Approx) intersects the boundary and use this
-- point as (x2Approx,y2Approx).
-- If the resulting point is on the left or right boundary,
-- we will now consider x as the 'incremented variable' and we
-- will compute the y-coordinate using a Newton iteration.
-- Similarly, if the point is on the top or bottom boundary,
-- we will consider y as the 'incremented variable' and we
-- will compute the x-coordinate using a Newton iteration.
if x2Approx >= xMaxSF then
incVar := x
lookingFor := BDRY
x2Approx := xMaxSF
y2Approx := y1 + (-px/py)*(x2Approx - x1)
else
if x2Approx <= xMinSF then
incVar := x
lookingFor := BDRY
x2Approx := xMinSF
y2Approx := y1 + (-px/py)*(x2Approx - x1)
if y2Approx >= yMaxSF then
incVar := y
lookingFor := BDRY
y2Approx := yMaxSF
x2Approx := x1 + (-py/px)*(y2Approx - y1)
else
if y2Approx <= yMinSF then
incVar := y
lookingFor := BDRY
y2Approx := yMinSF
x2Approx := x1 + (-py/px)*(y2Approx - y1)
-- set xLo = min(x1,x2Approx), xHi = max(x1,x2Approx)
-- set yLo = min(y1,y2Approx), yHi = max(y1,y2Approx)
if x1 < x2Approx then
xLo := x1
xHi := x2Approx
else
xLo := x2Approx
xHi := x1
if y1 < y2Approx then
yLo := y1
yHi := y2Approx
else
yLo := y2Approx
yHi := y1
-- check for critical points (x*,y*) with x* between
-- x1 and x2Approx or y* between y1 and y2Approx
-- store values of x2Approx and y2Approx
x2Approxx := x2Approx
y2Approxx := y2Approx
-- xPointList will contain all critical points (x*,y*)
-- with x* between x1 and x2Approx
xPointList : L Pt := nil()
-- yPointList will contain all critical points (x*,y*)
-- with y* between y1 and y2Approx
yPointList : L Pt := nil()
for pt in crits repeat
xx := xCoord pt; yy := yCoord pt
-- if x1 = x2Approx, then p1 is a point with horizontal
-- tangent line
-- in this case, we don't want critical points with
-- x-coordinate x1
if xx = x2Approx and not (xx = x1) then
if min(abs(yy-yLo),abs(yy-yHi)) < delta then
xPointList := cons(pt,xPointList)
if ((xLo < xx) and (xx < xHi)) then
if min(abs(yy-yLo),abs(yy-yHi)) < delta then
xPointList := cons(pt,nil())
x2Approx := xx
if xx < x1 then xLo := xx else xHi := xx
-- if y1 = y2Approx, then p1 is a point with vertical
-- tangent line
-- in this case, we don't want critical points with
-- y-coordinate y1
if yy = y2Approx and not (yy = y1) then
yPointList := cons(pt,yPointList)
if ((yLo < yy) and (yy < yHi)) then
if min(abs(xx-xLo),abs(xx-xHi)) < delta then
yPointList := cons(pt,nil())
y2Approx := yy
if yy < y1 then yLo := yy else yHi := yy
-- points in both xPointList and yPointList
if (not null xPointList) and (not null yPointList) then
xPointList = yPointList =>
-- this implies that the lists have only one point
incVar := incVar0
if incVar = x then
y2Approx := y1 + (-px/py)*(x2Approx - x1)
else
x2Approx := x1 + (-py/px)*(y2Approx - y1)
lookingFor := CRIT -- proceed
incVar0 = x =>
-- first try Newton iteration with 'y' as incremented variable
x2Temp := x1 + (-py/px)*(y2Approx - y1)
f := SFPolyToUPoly(eval(pSF,y,y2Approx))
x2New := newtonApprox(f,x2Temp,err,bound)
x2New case "failed" =>
y2Approx := y1 + (-px/py)*(x2Approx - x1)
incVar := x
lookingFor := CRIT -- proceed
y2Temp := y1 + (-px/py)*(x2Approx - x1)
f := SFPolyToUPoly(eval(pSF,x,x2Approx))
y2New := newtonApprox(f,y2Temp,err,bound)
y2New case "failed" =>
return computeNextPt(pSF,dpdxSF,dpdySF,x,y,p0,p1,corners,_
abs((x2Approx-x1)/2),err,bound,crits,bdry)
pt1 := makePt(x2Approx,y2New :: SF)
pt2 := makePt(x2New :: SF,y2Approx)
critPt1 := findPtOnList(pt1,crits)
critPt2 := findPtOnList(pt2,crits)
(critPt1 case "failed") and (critPt2 case "failed") =>
abs(x2Approx - x1) > abs(x2Temp - x1) =>
return [pt1,NADA]
return [pt2,NADA]
(critPt1 case "failed") =>
return [critPt2::Pt,CRIT]
(critPt2 case "failed") =>
return [critPt1::Pt,CRIT]
abs(x2Approx - x1) > abs(x2Temp - x1) =>
return [critPt2::Pt,CRIT]
return [critPt1::Pt,CRIT]
y2Temp := y1 + (-px/py)*(x2Approx - x1)
f := SFPolyToUPoly(eval(pSF,x,x2Approx))
y2New := newtonApprox(f,y2Temp,err,bound)
y2New case "failed" =>
x2Approx := x1 + (-py/px)*(y2Approx - y1)
incVar := y
lookingFor := CRIT -- proceed
x2Temp := x1 + (-py/px)*(y2Approx - y1)
f := SFPolyToUPoly(eval(pSF,y,y2Approx))
x2New := newtonApprox(f,x2Temp,err,bound)
x2New case "failed" =>
return computeNextPt(pSF,dpdxSF,dpdySF,x,y,p0,p1,corners,_
abs((y2Approx-y1)/2),err,bound,crits,bdry)
pt1 := makePt(x2Approx,y2New :: SF)
pt2 := makePt(x2New :: SF,y2Approx)
critPt1 := findPtOnList(pt1,crits)
critPt2 := findPtOnList(pt2,crits)
(critPt1 case "failed") and (critPt2 case "failed") =>
abs(y2Approx - y1) > abs(y2Temp - y1) =>
return [pt2,NADA]
return [pt1,NADA]
(critPt1 case "failed") =>
return [critPt2::Pt,CRIT]
(critPt2 case "failed") =>
return [critPt1::Pt,CRIT]
abs(y2Approx - y1) > abs(y2Temp - y1) =>
return [critPt1::Pt,CRIT]
return [critPt2::Pt,CRIT]
if (not null xPointList) and (null yPointList) then
y2Approx := y1 + (-px/py)*(x2Approx - x1)
incVar0 = x =>
incVar := x
lookingFor := CRIT -- proceed
f := SFPolyToUPoly(eval(pSF,x,x2Approx))
y2New := newtonApprox(f,y2Approx,err,bound)
y2New case "failed" =>
x2Approx := x2Approxx
y2Approx := y2Approxx -- proceed
pt := makePt(x2Approx,y2New::SF)
critPt := findPtOnList(pt,crits)
critPt case "failed" =>
return [pt,NADA]
return [critPt :: Pt,CRIT]
if (null xPointList) and (not null yPointList) then
x2Approx := x1 + (-py/px)*(y2Approx - y1)
incVar0 = y =>
incVar := y
lookingFor := CRIT -- proceed
f := SFPolyToUPoly(eval(pSF,y,y2Approx))
x2New := newtonApprox(f,x2Approx,err,bound)
x2New case "failed" =>
x2Approx := x2Approxx
y2Approx := y2Approxx -- proceed
pt := makePt(x2New::SF,y2Approx)
critPt := findPtOnList(pt,crits)
critPt case "failed" =>
return [pt,NADA]
return [critPt :: Pt,CRIT]
if incVar = x then
x2 := x2Approx
f := SFPolyToUPoly(eval(pSF,x,x2))
y2New := newtonApprox(f,y2Approx,err,bound)
y2New case "failed" =>
return computeNextPt(pSF,dpdxSF,dpdySF,x,y,p0,p1,corners,_
abs((x2-x1)/2),err,bound,crits,bdry)
y2 := y2New :: SF
else
y2 := y2Approx
f := SFPolyToUPoly(eval(pSF,y,y2))
x2New := newtonApprox(f,x2Approx,err,bound)
x2New case "failed" =>
return computeNextPt(pSF,dpdxSF,dpdySF,x,y,p0,p1,corners,_
abs((y2-y1)/2),err,bound,crits,bdry)
x2 := x2New :: SF
pt := makePt(x2,y2)
--!! check that 'pt' is not out of bounds
-- check if you've gotten a critical or boundary point
lookingFor = NADA =>
[pt,lookingFor]
lookingFor = BDRY =>
bdryPt := findPtOnList(pt,bdry)
bdryPt case "failed" =>
error "couldn't find boundary point"
[bdryPt :: Pt,BDRY]
critPt := findPtOnList(pt,crits)
critPt case "failed" =>
[pt,NADA]
[critPt :: Pt,CRIT]
--% Newton iterations
newtonApprox(f,a0,err,bound) ==
-- Newton iteration to approximate a root of the polynomial 'f'
-- using an initial approximation of 'a0'
-- Newton iteration terminates when consecutive approximations
-- are within 'err' of each other
-- returns "failed" if this has not been achieved after 'bound'
-- iterations
Df := differentiate f
oldApprox := a0
newApprox := a0 - elt(f,a0)/elt(Df,a0)
i : PI := 1
while abs(newApprox - oldApprox) > err repeat
i = bound => return "failed"
oldApprox := newApprox
newApprox := oldApprox - elt(f,oldApprox)/elt(Df,oldApprox)
i := i+1
newApprox
--% graphics output
listBranches(acplot) == acplot.branches
--% terminal output
coerce(acplot:%) ==
pp := acplot.poly :: OUT
xx := acplot.xVar :: OUT
yy := acplot.yVar :: OUT
xLo := acplot.minXVal :: OUT
xHi := acplot.maxXVal :: OUT
yLo := acplot.minYVal :: OUT
yHi := acplot.maxYVal :: OUT
zip := message(" = 0")
com := message(", ")
les := message(" <= ")
l : L OUT :=
[pp,zip,com,xLo,les,xx,les,xHi,com,yLo,les,yy,les,yHi]
f : L OUT := nil()
for branch in acplot.branches repeat
ll : L OUT := [p :: OUT for p in branch]
f := cons(vconcat ll,f)
ff := vconcat(hconcat l,vconcat f)
vconcat(message "ACPLOT",ff)
@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
-- - Redistributions of source code must retain the above copyright
-- notice, this list of conditions and the following disclaimer.
--
-- - Redistributions in binary form must reproduce the above copyright
-- notice, this list of conditions and the following disclaimer in
-- the documentation and/or other materials provided with the
-- distribution.
--
-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the
-- names of its contributors may be used to endorse or promote products
-- derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>
<<package REALSOLV RealSolvePackage>>
<<domain ACPLOT PlaneAlgebraicCurvePlot>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}