\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra realzero.spad}
\author{Andy Neff}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package REAL0 RealZeroPackage}
<<package REAL0 RealZeroPackage>>=
)abbrev package REAL0 RealZeroPackage
++ Author: Andy Neff
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors: UnivariatePolynomial, RealZeroPackageQ
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This package provides functions for finding the real zeros
++ of univariate polynomials over the integers to arbitrary user-specified
++ precision. The results are returned as a list of
++ isolating intervals which are expressed as records with "left" and "right" rational number
++ components.

RealZeroPackage(Pol): T == C where
   Pol: UnivariatePolynomialCategory Integer
   RN ==> Fraction Integer
   Interval ==> Record(left : RN, right : RN)
   isoList ==> List(Interval)
   T == with
    -- next two functions find isolating intervals
      realZeros: (Pol) -> isoList
        ++ realZeros(pol) returns a list of isolating intervals for
        ++ all the real zeros of the univariate polynomial pol.
      realZeros: (Pol, Interval) -> isoList
        ++ realZeros(pol, range) returns a list of isolating intervals
        ++ for all the real zeros of the univariate polynomial pol which
        ++ lie in the interval expressed by the record range.
    -- next two functions return intervals smaller then tolerence
      realZeros: (Pol, RN) -> isoList
        ++ realZeros(pol, eps) returns a list of intervals of length less
        ++ than the rational number eps for all the real roots of the
        ++ polynomial pol.
      realZeros: (Pol, Interval, RN) -> isoList
        ++ realZeros(pol, int, eps) returns a list of intervals of length
        ++ less than the rational number eps for all the real roots of the
        ++ polynomial pol which lie in the interval expressed by the
        ++ record int.
      refine: (Pol, Interval, RN) -> Interval
        ++ refine(pol, int, eps) refines the interval int containing
        ++ exactly one root of the univariate polynomial pol to size less
        ++ than the rational number eps.
      refine: (Pol, Interval, Interval) -> Union(Interval,"failed")
        ++ refine(pol, int, range) takes a univariate polynomial pol and
        ++ and isolating interval int containing exactly one real
        ++ root of pol; the operation returns an isolating interval which
        ++ is contained within range, or "failed" if no such isolating interval exists.
      midpoint: Interval -> RN
        ++ midpoint(int) returns the midpoint of the interval int.
      midpoints: isoList -> List RN
        ++ midpoints(isolist) returns the list of midpoints for the list
        ++ of intervals isolist.
   C == add
      --Local Functions
      makeSqfr: Pol -> Pol
      ReZeroSqfr: (Pol) -> isoList
      PosZero: (Pol) -> isoList
      Zero1: (Pol) -> isoList
      transMult: (Integer, Pol) -> Pol
      transMultInv: (Integer, Pol) -> Pol
      transAdd1: (Pol) -> Pol
      invert: (Pol) -> Pol
      minus: (Pol) -> Pol
      negate: Interval -> Interval
      rootBound: (Pol) -> Integer
      var: (Pol) -> Integer

      negate(int : Interval):Interval == [-int.right,-int.left]

      midpoint(i : Interval):RN ==  (1/2)*(i.left + i.right)

      midpoints(li : isoList) : List RN ==
        [midpoint x for x in li]

      makeSqfr(F : Pol):Pol ==
         sqfr := squareFree F
         F := */[s.factor for s in factors(sqfr)]

      realZeros(F : Pol) ==
         ReZeroSqfr makeSqfr F

      realZeros(F : Pol, rn : RN) ==
         F := makeSqfr F
         [refine(F,int,rn) for int in ReZeroSqfr(F)]

      realZeros(F : Pol, bounds : Interval) ==
         F := makeSqfr F
         [rint::Interval for int in ReZeroSqfr(F) |
             (rint:=refine(F,int,bounds)) case Interval]

      realZeros(F : Pol, bounds : Interval, rn : RN) ==
         F := makeSqfr F
         [refine(F,int,rn) for int in realZeros(F,bounds)]

      ReZeroSqfr(F : Pol) ==
         F = 0 => error "ReZeroSqfr: zero polynomial"
         L : isoList := []
         degree(F) = 0 => L
         if positive?(r := minimumDegree(F)) then
              L := [[0,0]$Interval]
              tempF := F exquo monomial(1, r)
              if not (tempF case "failed") then
                   F := tempF
         J:isoList := [negate int for int in reverse(PosZero(minus(F)))]
         K : isoList := PosZero(F)
         append(append(J, L), K)

      PosZero(F : Pol) ==   --F is square free, primitive
                          --and F(0) ~= 0; returns isoList for positive
                          --roots of F

         b : Integer := rootBound(F)
         F := transMult(b,F)
         L : isoList := Zero1(F)
         int : Interval
         L := [[b*int.left, b*int.right]$Interval for int in L]

      Zero1(F : Pol) ==   --returns isoList for roots of F in (0,1)
         J : isoList
         K : isoList
         L : isoList
         L := []
         (v := var(transAdd1(invert(F)))) = 0 => []
         v = 1 => L := [[0,1]$Interval]
         G : Pol := transMultInv(2, F)
         H : Pol := transAdd1(G)
         if positive? minimumDegree H then
                 -- H has a root at 0 => F has one at 1/2, and G at 1
              L := [[1/2,1/2]$Interval]
              Q : Pol := monomial(1, 1)
              tempH : Union(Pol, "failed") := H exquo Q
              if not (tempH case "failed") then H := tempH
              Q := Q + monomial(-1, 0)
              tempG : Union(Pol, "failed") := G exquo Q
              if not (tempG case "failed") then G := tempG
         int : Interval
         J := [[(int.left+1)* (1/2),(int.right+1) * (1/2)]$Interval
                                                 for int in Zero1(H)]
         K := [[int.left * (1/2), int.right * (1/2)]$Interval
                                                 for int in Zero1(G)]
         append(append(J, L), K)

      rootBound(F : Pol) ==  --returns power of 2 that is a bound
                             --for the positive roots of F
         if negative? leadingCoefficient(F) then F := -F
         lcoef := leadingCoefficient(F)
         F := reductum(F)
         i : Integer := 0
         while not (F = 0) repeat
              if negative?(an := leadingCoefficient(F)) then i := i - an
              F := reductum(F)
         b : Integer := 1
         while (b * lcoef) <= i repeat
              b := 2 * b
         b

      transMult(c : Integer, F : Pol) ==
                  --computes Pol G such that G(x) = F(c*x)
         G : Pol := 0
         while not (F = 0) repeat
              n := degree(F)
              G := G + monomial((c**n) * leadingCoefficient(F), n)
              F := reductum(F)
         G

      transMultInv(c : Integer, F : Pol) ==
                  --computes Pol G such that G(x) = (c**n) * F(x/c)
         d := degree(F)
         cc : Integer := 1
         G : Pol := monomial(leadingCoefficient F,d)
         while (F:=reductum(F)) ~= 0 repeat
              n := degree(F)
              cc := cc*(c**(d-n):NonNegativeInteger)
              G := G + monomial(cc * leadingCoefficient(F), n)
              d := n
         G

--    otransAdd1(F : Pol) ==
--                --computes Pol G such that G(x) = F(x+1)
--       G : Pol := F
--       n : Integer := 1
--       while (F := differentiate(F)) ~= 0 repeat
--            if not ((tempF := F exquo n) case "failed") then F := tempF
--            G := G + F
--            n := n + 1
--       G

      transAdd1(F : Pol) ==
                  --computes Pol G such that G(x) = F(x+1)
         n := degree F
         v := vectorise(F, n+1)
         for i in 0..(n-1) repeat
            for j in (n-i)..n repeat
               qsetelt!(v,j, qelt(v,j) + qelt(v,(j+1)))
         ans : Pol := 0
         for i in 0..n repeat
            ans := ans + monomial(qelt(v,(i+1)),i)
         ans


      minus(F : Pol) ==
                  --computes Pol G such that G(x) = F(-x)
         G : Pol := 0
         while not (F = 0) repeat
              n := degree(F)
              coef := leadingCoefficient(F)
              odd? n =>
                   G := G + monomial(-coef, n)
                   F := reductum(F)
              G := G + monomial(coef, n)
              F := reductum(F)
         G

      invert(F : Pol) ==
                  --computes Pol G such that G(x) = (x**n) * F(1/x)
         G : Pol := 0
         n := degree(F)
         while not (F = 0) repeat
              G := G + monomial(leadingCoefficient(F),
                                (n-degree(F))::NonNegativeInteger)
              F := reductum(F)
         G

      var(F : Pol) ==    --number of sign variations in coefs of F
         i : Integer := 0
         LastCoef : Boolean
         next : Boolean
         LastCoef := negative? leadingCoefficient(F)
         while not ((F := reductum(F)) = 0) repeat
              next := negative? leadingCoefficient(F)
              if ((not LastCoef) and next) or
                  ((not next) and LastCoef) then i := i+1
              LastCoef := next
         i

      refine(F : Pol, int : Interval, bounds : Interval) ==
         lseg := min(int.right,bounds.right) - max(int.left,bounds.left)
         negative? lseg => "failed"
         lseg = 0 =>
            pt :=
               int.left = bounds.right => int.left
               int.right
            elt(transMultInv(denom(pt),F),numer pt) = 0 => [pt,pt]
            "failed"
         lseg = int.right - int.left => int
         refine(F, refine(F, int, lseg), bounds)

      refine(F : Pol, int : Interval, eps : RN) ==
         a := int.left
         b := int.right
         a=b => [a,b]$Interval
         an : Integer := numer(a)
         ad : Integer := denom(a)
         bn : Integer := numer(b)
         bd : Integer := denom(b)
         xfl : Boolean := false
         if (u:=elt(transMultInv(ad, F), an)) = 0 then
             F := (F exquo (monomial(ad,1)-monomial(an,0)))::Pol
             u:=elt(transMultInv(ad, F), an)
         if (v:=elt(transMultInv(bd, F), bn)) = 0 then
             F := (F exquo (monomial(bd,1)-monomial(bn,0)))::Pol
             v:=elt(transMultInv(bd, F), bn)
             u:=elt(transMultInv(ad, F), an)
         if positive? u then (F:=-F;v:=-v)
         if negative? v then
            error [int, "is not a valid isolation interval for", F]
         if eps <= 0 then error "precision must be positive"
         while (b - a) >= eps repeat
              mid : RN := (b + a) * (1/2)
              midn : Integer := numer(mid)
              midd : Integer := denom(mid)
              negative?(v := elt(transMultInv(midd, F), midn)) =>
                   a := mid
                   an := midn
                   ad := midd
              positive? v =>
                   b := mid
                   bn := midn
                   bd := midd
              v = 0 =>
                   a := mid
                   b := mid
                   an := midn
                   ad := midd
                   xfl := true
         [a, b]$Interval

@
\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 REAL0 RealZeroPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}