axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / ACPLOT.spad

)abbrev domain ACPLOT PlaneAlgebraicCurvePlot
++ Author: Clifton J. Williamson and Timothy Daly
++ Date Created: Fall 1988
++ Date Last Updated: 27 April 1990
++ Description:
++ Plot a NON-SINGULAR plane algebraic curve p(x,y) = 0.

PlaneAlgebraicCurvePlot() : SIG == CODE where

  SIG ==> PlottablePlaneCurveCategory with

    makeSketch : (Polynomial Integer,Symbol,Symbol,Segment Fraction Integer,
                 Segment Fraction Integer) -> %
      ++ makeSketch(p,x,y,a..b,c..d) creates an ACPLOT of the
      ++ curve \spad{p = 0} in the region a <= x <= b, c <= y <= d.
      ++ More specifically, 'makeSketch' plots a non-singular algebraic curve
      ++ \spad{p = 0} in an rectangular region xMin <= x <= xMax,
      ++ 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.
      ++
      ++X makeSketch(x+y,x,y,-1/2..1/2,-1/2..1/2)$ACPLOT

    refine : (%,DoubleFloat) -> %
      ++ refine(p,x) is not documented
      ++
      ++X sketch:=makeSketch(x+y,x,y,-1/2..1/2,-1/2..1/2)$ACPLOT
      ++X refined:=refine(sketch,0.1)

  CODE ==> add

    import PointPackage DoubleFloat
    import Plot
    import RealSolvePackage

    BoundaryPts ==> Record(left:   List Point DoubleFloat,_
                           right:  List Point DoubleFloat,_
                           bottom: List Point DoubleFloat,_
                           top:    List Point DoubleFloat)

    NewPtInfo   ==> Record(newPt: Point DoubleFloat,_
                           type:  String)

    Corners     ==> Record(minXVal: DoubleFloat,_
                           maxXVal: DoubleFloat,_
                           minYVal: DoubleFloat,_
                           maxYVal: DoubleFloat)

    kinte       ==> solve$RealSolvePackage()
  
    rsolve      ==> realSolve$RealSolvePackage()
  
    singValBetween?:(DoubleFloat,DoubleFloat,List DoubleFloat) -> Boolean
  
    segmentInfo:(DoubleFloat -> DoubleFloat,DoubleFloat,DoubleFloat,_
                 List DoubleFloat,List DoubleFloat,List DoubleFloat,_
                 DoubleFloat,DoubleFloat) -> _
      Record(seg:Segment DoubleFloat,_
             left: DoubleFloat,_
             lowerVals: List DoubleFloat,_
             upperVals:List DoubleFloat)
  
    swapCoords:Point DoubleFloat -> Point DoubleFloat
  
    samePlottedPt?:(Point DoubleFloat,Point DoubleFloat) -> Boolean
  
    findPtOnList:(Point DoubleFloat,List Point DoubleFloat) -> _
      Union(Point DoubleFloat,"failed")
  
    makeCorners:(DoubleFloat,DoubleFloat,DoubleFloat,DoubleFloat) -> Corners
  
    getXMin: Corners -> DoubleFloat
  
    getXMax: Corners -> DoubleFloat
  
    getYMin: Corners -> DoubleFloat
  
    getYMax: Corners -> DoubleFloat
  
    SFPolyToUPoly:Polynomial DoubleFloat -> _
      SparseUnivariatePolynomial DoubleFloat
  
    RNPolyToUPoly:Polynomial Fraction Integer -> _
      SparseUnivariatePolynomial Fraction Integer
  
    coerceCoefsToSFs:Polynomial Integer -> Polynomial DoubleFloat
  
    coerceCoefsToRNs:Polynomial Integer -> Polynomial Fraction Integer
  
    RNtoSF:Fraction Integer -> DoubleFloat
  
    RNtoNF:Fraction Integer -> Float
  
    SFtoNF:DoubleFloat -> Float
  
    listPtsOnHorizBdry:(Polynomial Fraction Integer,Symbol,Fraction Integer,_
                        Float,Float) -> _
      List Point DoubleFloat
  
    listPtsOnVertBdry:(Polynomial Fraction Integer,Symbol,Fraction Integer,_
                       Float,Float) -> _
      List Point DoubleFloat
  
    listPtsInRect:(List List Float,Float,Float,Float,Float) -> _
      List Point DoubleFloat
  
    ptsSuchThat?:(List List Float,List Float -> Boolean) -> Boolean
  
    inRect?:(List Float,Float,Float,Float,Float) -> Boolean
  
    onHorzSeg?:(List Float,Float,Float,Float) -> Boolean
  
    onVertSeg?:(List Float,Float,Float,Float) -> Boolean
  
    newX:(List List Float,List List Float,Float,Float,Float,Fraction Integer,_
          Fraction Integer) -> Fraction Integer
  
    newY:(List List Float,List List Float,Float,Float,Float,_
          Fraction Integer,Fraction Integer) -> Fraction Integer
  
    makeOneVarSketch:(Polynomial Integer,Symbol,Symbol,Fraction Integer,_
                      Fraction Integer,Fraction Integer,Fraction Integer,_
                      Symbol) -> %
  
    makeLineSketch:(Polynomial Integer,Symbol,Symbol,Fraction Integer,_
                    Fraction Integer,Fraction Integer,Fraction Integer) -> %
  
    makeRatFcnSketch:(Polynomial Integer,Symbol,Symbol,Fraction Integer,_
                      Fraction Integer,Fraction Integer,Fraction Integer,_
                      Symbol) -> %
  
    makeGeneralSketch:(Polynomial Integer,Symbol,Symbol,Fraction Integer,_
                       Fraction Integer,Fraction Integer,Fraction Integer) -> %
  
    traceBranches:(Polynomial DoubleFloat,Polynomial DoubleFloat,_
                   Polynomial DoubleFloat,Symbol,Symbol,Corners,DoubleFloat,_
                   DoubleFloat,PositiveInteger, List Point DoubleFloat,_
                   BoundaryPts) -> List List Point DoubleFloat
  
    dummyFirstPt:(Point DoubleFloat,Polynomial DoubleFloat,_
                  Polynomial DoubleFloat,Symbol,Symbol,List Point DoubleFloat,_
                  List Point DoubleFloat,List Point DoubleFloat,_
                  List Point DoubleFloat) -> Point DoubleFloat
  
    listPtsOnSegment:(Polynomial DoubleFloat,Polynomial DoubleFloat,_
                      Polynomial DoubleFloat,Symbol,Symbol,Point DoubleFloat,_
                      Point DoubleFloat,Corners, DoubleFloat,DoubleFloat,_
                      PositiveInteger,List Point DoubleFloat,_
                      List Point DoubleFloat) -> List List Point DoubleFloat
  
    listPtsOnLoop:(Polynomial DoubleFloat,Polynomial DoubleFloat,_
                   Polynomial DoubleFloat,Symbol,Symbol,Point DoubleFloat,_
                   Corners, DoubleFloat,DoubleFloat,PositiveInteger,_
                   List Point DoubleFloat,List Point DoubleFloat) -> _
                   List List Point DoubleFloat
  
    computeNextPt:(Polynomial DoubleFloat,Polynomial DoubleFloat,_
                   Polynomial DoubleFloat,Symbol,Symbol,Point DoubleFloat,_
                   Point DoubleFloat,Corners, DoubleFloat,DoubleFloat,_
                   PositiveInteger,List Point DoubleFloat,_
                   List Point DoubleFloat) -> NewPtInfo
  
    newtonApprox:(SparseUnivariatePolynomial DoubleFloat, DoubleFloat, _
                  DoubleFloat, PositiveInteger) -> Union(DoubleFloat, "failed")
  
  --% representation
  
    Rep := Record(poly    : Polynomial Integer,_
                  xVar    : Symbol,_
                  yVar    : Symbol,_
                  minXVal : Fraction Integer,_
                  maxXVal : Fraction Integer,_
                  minYVal : Fraction Integer,_
                  maxYVal : Fraction Integer,_
                  bdryPts : BoundaryPts,_
                  hTanPts : List Point DoubleFloat,_
                  vTanPts : List Point DoubleFloat,_
                  branches: List List Point DoubleFloat)
  
  --% global constants
  
    EPSILON : Float := .000001 -- precision to which realSolve finds roots
    PLOTERR : DoubleFloat := 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: Float -> DoubleFloat
    NFtoSF x == 0 + convert(x)$Float
  
  --% points
    makePt: (DoubleFloat,DoubleFloat) -> Point DoubleFloat
    makePt(xx,yy) == point(l : List DoubleFloat := [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(Integer,DoubleFloat)
  
    coerceCoefsToRNs(p) ==
    -- coefficients of 'p' are coerced to be DoubleFloat's
      map(coerce,p)$PolynomialFunctions2(Integer,Fraction Integer)
  
    RNtoSF(r) == coerce(r)@DoubleFloat
    RNtoNF(r) == coerce(r)@Float
    SFtoNF(x) == convert(x)@Float
  
  --% computation of special points
  
    listPtsOnHorizBdry(pRN,y,y0,xMinNF,xMaxNF) ==
    -- strict inequality here: corners on vertical boundary
      pointList : List Point DoubleFloat := nil()
      ySF := RNtoSF(y0)
      f := eval(pRN,y,y0)
      roots : List Float := 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 : List Point DoubleFloat := nil()
      xSF := RNtoSF(x0)
      f := eval(pRN,x,x0)
      roots : List Float := 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 : List Point DoubleFloat := 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 :List Float -> Boolean := x +-> inRect?(x,xLftNF,xRtNF,yMinNF,yMaxNF)
      ptsSuchThat?(singPts,foo) =>
        newX(vtanPts,singPts,yMinNF,yMaxNF,xNF,xRN,_
          horizInc/2::(Fraction Integer))
  --  ptsSuchThat?(vtanPts,onVertSeg?(#1,yMinNF,yMaxNF,xNewNF)) =>
      goo : List Float -> Boolean := x +-> onVertSeg?(x,yMinNF,yMaxNF,xNewNF)
      ptsSuchThat?(vtanPts,goo) =>
        newX(vtanPts,singPts,yMinNF,yMaxNF,xNF,xRN,_
          horizInc/2::(Fraction Integer))
      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:List Float -> Boolean := x +-> inRect?(x,xMinNF,xMaxNF,yBotNF,yTopNF)
      ptsSuchThat?(singPts,foo) =>
        newY(htanPts,singPts,xMinNF,xMaxNF,yNF,yRN,_
          vertInc/2::(Fraction Integer))
  --  ptsSuchThat?(htanPts,onHorzSeg?(#1,xMinNF,xMaxNF,yNewNF)) =>
      goo : List Float -> Boolean := x +-> onHorzSeg?(x,xMinNF,xMaxNF,yNewNF)
      ptsSuchThat?(htanPts,goo) =>
        newY(htanPts,singPts,xMinNF,xMaxNF,yNF,yRN,_
          vertInc/2::(Fraction Integer))
      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 : List Point DoubleFloat := nil()
      rt : List Point DoubleFloat := nil()
      bt : List Point DoubleFloat := nil() 
      tp : List Point DoubleFloat := nil()
      htans : List Point DoubleFloat := nil() 
      vtans : List Point DoubleFloat := nil()
      bran : List List Point DoubleFloat := nil()
      roots := kinte(p,EPSILON)
      sketchRoots : List DoubleFloat := 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 : List Point DoubleFloat := [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 : List Point DoubleFloat := [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)) :: DoubleFloat
      b := ground(coefficient(p,y,1)) :: DoubleFloat
      c := ground(coefficient(coefficient(p,x,0),y,0)) :: DoubleFloat
      lf : List Point DoubleFloat := nil()
      rt : List Point DoubleFloat := nil()
      bt : List Point DoubleFloat := nil()
      tp : List Point DoubleFloat := nil()
      htans : List Point DoubleFloat := nil()
      vtans : List Point DoubleFloat := nil()
      branch : List Point DoubleFloat := nil()
      bran : List List Point DoubleFloat := 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::DoubleFloat)
            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: for example, upper left corner is on graph with y' > 0
      factoredP := p ::(Factored Polynomial Integer)
      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:DoubleFloat -> DoubleFloat := s +-> elt(numUPolySF,s)/elt(denUPolySF,s)
      -- 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 : List Point DoubleFloat := nil()
      rt : List Point DoubleFloat := nil()
      bt : List Point DoubleFloat := nil()
      tp : List Point DoubleFloat := 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 : List DoubleFloat :=
        roots : List Float := kinte(den,EPSILON)
        outList : List DoubleFloat := nil()
        for root in roots repeat
          if (indVarMinNF < root) and (root < indVarMaxNF) then
            outList := cons(NFtoSF root,outList)
        sort((x,y) +-> x < y, outList)
      topList : List DoubleFloat :=
        roots : List Float := kinte(eval(pRN,depVar,depVarMax),EPSILON)
        outList : List DoubleFloat := nil()
        for root in roots repeat
          if (indVarMinNF < root) and (root < indVarMaxNF) then
            outList := cons(NFtoSF root,outList)
        sort((x,y) +-> x < y, outList)
      botList : List DoubleFloat :=
        roots : List Float := kinte(eval(pRN,depVar,depVarMin),EPSILON)
        outList : List DoubleFloat := nil()
        for root in roots repeat
          if (indVarMinNF < root) and (root < indVarMaxNF) then
            outList := cons(NFtoSF root,outList)
        sort((x,y) +-> x < y, 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$(Fraction Integer) then
        leftVal  := depVarMinSF - (abs(depVarMinSF) + 1$DoubleFloat)
      else
        leftVal  := f(indVarMinSF)
      if elt(denUPolyRN,indVarMax) = 0$(Fraction Integer) then
        rightVal := depVarMinSF - (abs(depVarMinSF) + 1$DoubleFloat)
      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 : List List Point DoubleFloat := 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$DoubleFloat  -- @#$%^&* 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 : List Segment DoubleFloat := 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 :: (Factored Polynomial Integer)
      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:List Float -> Boolean := s +-> inRect?(s,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 : Fraction Integer := (1/20) * (xMax - xMin)
      vInc : Fraction Integer := (1/20) * (yMax - yMin)
  --  if ptsSuchThat?(critPts,onVertSeg?(#1,yMinNF,yMaxNF,xMinNF)) then
      foo : List Float -> Boolean := s +-> onVertSeg?(s,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 : List Float -> Boolean := s +-> onVertSeg?(s,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 : List Float -> Boolean := s +-> onHorzSeg?(s,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 : List Float -> Boolean := s +-> onHorzSeg?(s,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 : PositiveInteger := 10
      -- compute points on the boundary
      pRN := coerceCoefsToRNs(p)
      lf : List Point DoubleFloat := 
        listPtsOnVertBdry(pRN,x,xMin,yMinNF,yMaxNF)
      rt : List Point DoubleFloat := 
        listPtsOnVertBdry(pRN,x,xMax,yMinNF,yMaxNF)
      bt : List Point DoubleFloat := 
        listPtsOnHorizBdry(pRN,y,yMin,xMinNF,xMaxNF)
      tp : List Point DoubleFloat := 
        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 : PositiveInteger := 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 : List List Point DoubleFloat := 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$DoubleFloat; one := 1$DoubleFloat
      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 : List Point DoubleFloat := [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$DoubleFloat,y1 - 1$DoubleFloat)
      pointList : List Point DoubleFloat := [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:Point DoubleFloat,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 : List Point DoubleFloat := nil()
      -- yPointList will contain all critical points (x*,y*)
      -- with y* between y1 and y2Approx
      yPointList : List Point DoubleFloat := 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 :: DoubleFloat)
          pt2 := makePt(x2New :: DoubleFloat,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::(Point DoubleFloat),CRIT]
          (critPt2 case "failed") =>
            return [critPt1::(Point DoubleFloat),CRIT]
          abs(x2Approx - x1) > abs(x2Temp - x1) =>
            return [critPt2::(Point DoubleFloat),CRIT]
          return [critPt1::(Point DoubleFloat),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 :: DoubleFloat)
        pt2 := makePt(x2New :: DoubleFloat,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::(Point DoubleFloat),CRIT]
        (critPt2 case "failed") =>
          return [critPt1::(Point DoubleFloat),CRIT]
        abs(y2Approx - y1) > abs(y2Temp - y1) =>
          return [critPt1::(Point DoubleFloat),CRIT]
        return [critPt2::(Point DoubleFloat),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::DoubleFloat)
        critPt := findPtOnList(pt,crits)
        critPt case "failed" =>
          return [pt,NADA]
        return [critPt :: (Point DoubleFloat),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::DoubleFloat,y2Approx)
        critPt := findPtOnList(pt,crits)
        critPt case "failed" =>
          return [pt,NADA]
        return [critPt :: (Point DoubleFloat),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 :: DoubleFloat
      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 :: DoubleFloat
      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 :: (Point DoubleFloat),BDRY]
      critPt := findPtOnList(pt,crits)
      critPt case "failed" =>
        [pt,NADA]
      [critPt :: (Point DoubleFloat),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 : PositiveInteger := 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 :: OutputForm
      xx := acplot.xVar :: OutputForm
      yy := acplot.yVar :: OutputForm
      xLo := acplot.minXVal :: OutputForm
      xHi := acplot.maxXVal :: OutputForm
      yLo := acplot.minYVal :: OutputForm
      yHi := acplot.maxYVal :: OutputForm
      zip := message(" = 0")
      com := message(",   ")
      les := message(" <= ")
      l : List OutputForm :=
        [pp,zip,com,xLo,les,xx,les,xHi,com,yLo,les,yy,les,yHi]
      f : List OutputForm := nil()
      for branch in acplot.branches repeat
        ll : List OutputForm := [p :: OutputForm for p in branch]
        f := cons(vconcat ll,f)
      ff := vconcat(hconcat l,vconcat f)
      vconcat(message "ACPLOT",ff)