\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra sturm.spad}
\author{Lalo Gonzalez-Vega}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package SHP SturmHabichtPackage}
<<package SHP SturmHabichtPackage>>=
)abbrev package SHP SturmHabichtPackage
++ Author: Lalo Gonzalez-Vega
++ Date Created: 1994?
++ Date Last Updated: 30 January 1996
++ Basic Functions: 
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: localization
++ References:
++ Description: 
++ This package produces functions for counting
++ etc. real roots of univariate polynomials in x over R, which must
++ be an OrderedIntegralDomain
SturmHabichtPackage(R,x): T == C where
  R: OrderedIntegralDomain
  x: Symbol

  UP ==> UnivariatePolynomial
  L ==> List
  INT ==> Integer
  NNI ==> NonNegativeInteger

  T == with
--     subresultantSequenceBegin:(UP(x,R),UP(x,R)) -> L UP(x,R)
--       ++ \spad{subresultantSequenceBegin(p1,p2)} computes the  initial terms
--       ++ of the Subresultant sequence Sres(j)(P,deg(P),Q,deg(P)-1)
--       ++ when deg(Q)<deg(P)
--     subresultantSequenceNext:L UP(x,R) -> L UP(x,R)
--     subresultantSequenceInner:(UP(x,R),UP(x,R)) -> L UP(x,R)
     subresultantSequence:(UP(x,R),UP(x,R)) -> L UP(x,R)
       ++ subresultantSequence(p1,p2) computes the (standard) 
       ++ subresultant sequence of p1 and p2
--     sign:R -> R
--     delta:NNI -> R

--     polsth1:(UP(x,R),NNI,UP(x,R),NNI,R) -> L UP(x,R)
--     polsth2:(UP(x,R),NNI,UP(x,R),NNI,R) -> L UP(x,R)
--     polsth3:(UP(x,R),NNI,UP(x,R),NNI,R) -> L UP(x,R)
     SturmHabichtSequence:(UP(x,R),UP(x,R)) -> L UP(x,R)
       ++ SturmHabichtSequence(p1,p2) computes the Sturm-Habicht
       ++ sequence of p1 and p2
     SturmHabichtCoefficients:(UP(x,R),UP(x,R)) -> L R
       ++ SturmHabichtCoefficients(p1,p2) computes the principal
       ++ Sturm-Habicht coefficients of p1 and p2

--     variation:L R -> INT
--     permanence:L R -> INT
--     qzeros:L R -> L R
--     epsil:(NNI,R,R) -> INT
--     numbnce:L R -> NNI
--     numbce:L R -> NNI
--     wfunctaux:L R -> INT
--     wfunct:L R -> INT

     SturmHabicht:(UP(x,R),UP(x,R)) -> INT 
       ++ SturmHabicht(p1,p2) computes c_{+}-c_{-} where
       ++ c_{+} is the number of real roots of p1 with p2>0 and c_{-}
       ++ is the number of real roots of p1 with p2<0. If p2=1 what
       ++ you get is the number of real roots of p1.
     countRealRoots:(UP(x,R)) -> INT
       ++ countRealRoots(p) says how many real roots p has
     if R has GcdDomain then
        SturmHabichtMultiple:(UP(x,R),UP(x,R)) -> INT
          ++ SturmHabichtMultiple(p1,p2) computes c_{+}-c_{-} where
          ++ c_{+} is the number of real roots of p1 with p2>0 and c_{-}
          ++ is the number of real roots of p1 with p2<0. If p2=1 what
          ++ you get is the number of real roots of p1.
        countRealRootsMultiple:(UP(x,R)) -> INT
          ++ countRealRootsMultiple(p) says how many real roots p has,
          ++ counted with multiplicity


  C == add
     p1,p2: UP(x,R)
     Ex ==> OutputForm
     import OutputForm

     subresultantSequenceBegin(p1,p2):L UP(x,R) ==
       d1:NNI:=degree(p1)
       d2:NNI:=degree(p2)
       n:NNI:=(d1-1)::NNI
       d2 = n =>
         Pr:UP(x,R):=pseudoRemainder(p1,p2)
         append([p1,p2]::L UP(x,R),[Pr]::L UP(x,R))
       d2 = (n-1)::NNI =>
         Lc1:UP(x,R):=leadingCoefficient(p1)*leadingCoefficient(p2)*p2
         Lc2:UP(x,R):=-leadingCoefficient(p1)*pseudoRemainder(p1,p2)
         append([p1,p2]::L UP(x,R),[Lc1,Lc2]::L UP(x,R))
       LSubr:L UP(x,R):=[p1,p2]
       in1:INT:=(d2+1)::INT
       in2:INT:=(n-1)::INT
       for i in in1..in2 repeat
         LSubr:L UP(x,R):=append(LSubr::L UP(x,R),[0]::L UP(x,R))
       c1:R:=(leadingCoefficient(p1)*leadingCoefficient(p2))**((n-d2)::NNI)
       Lc1:UP(x,R):=monomial(c1,0)*p2
       Lc2:UP(x,R):=
         (-leadingCoefficient(p1))**((n-d2)::NNI)*pseudoRemainder(p1,p2)
       append(LSubr::L UP(x,R),[Lc1,Lc2]::L UP(x,R))

     subresultantSequenceNext(LcsI:L UP(x,R)):L UP(x,R) ==
       p2:UP(x,R):=last LcsI
       p1:UP(x,R):=first rest reverse LcsI
       d1:NNI:=degree(p1)
       d2:NNI:=degree(p2)
       in1:NNI:=(d1-1)::NNI
       d2 = in1 =>
         pr1:UP(x,R):=
           (pseudoRemainder(p1,p2) exquo (leadingCoefficient(p1))**2)::UP(x,R)
         append(LcsI:L UP(x,R),[pr1]:L UP(x,R))
       d2 < in1 =>
         c1:R:=leadingCoefficient(p1)
         pr1:UP(x,R):=
          (leadingCoefficient(p2)**((in1-d2)::NNI)*p2 exquo
              c1**((in1-d2)::NNI))::UP(x,R)
         pr2:UP(x,R):=
           (pseudoRemainder(p1,p2) exquo (-c1)**((in1-d2+2)::NNI))::UP(x,R)
         LSub:L UP(x,R):=[pr1,pr2]
         for k in ((d2+1)::INT)..((in1-1)::INT) repeat
           LSub:L UP(x,R):=append([0]:L UP(x,R),LSub:L UP(x,R))
         append(LcsI:L UP(x,R),LSub:L UP(x,R))

     subresultantSequenceInner(p1,p2):L UP(x,R) ==
       Lin:L UP(x,R):=subresultantSequenceBegin(p1:UP(x,R),p2:UP(x,R))
       indf:NNI:= if not(Lin.last::UP(x,R) = 0) then degree(Lin.last::UP(x,R))
                                               else 0
       while not(indf = 0) repeat
         Lin:L UP(x,R):=subresultantSequenceNext(Lin:L UP(x,R))
         indf:NNI:= if not(Lin.last::UP(x,R)=0) then degree(Lin.last::UP(x,R))
                                               else 0
       for j in #(Lin:L UP(x,R))..degree(p1) repeat
         Lin:L UP(x,R):=append(Lin:L UP(x,R),[0]:L UP(x,R))
       Lin


-- Computation of the subresultant sequence Sres(j)(P,p,Q,q) when:
--             deg(P) = p   and   deg(Q) = q   and   p > q

     subresultantSequence(p1,p2):L UP(x,R) ==
       p:NNI:=degree(p1)
       q:NNI:=degree(p2)
       List1:L UP(x,R):=subresultantSequenceInner(p1,p2)
       List2:L UP(x,R):=[p1,p2]
       c1:R:=leadingCoefficient(p1)
       for j in 3..#(List1) repeat
         Pr0:UP(x,R):=List1.j
         Pr1:UP(x,R):=(Pr0 exquo c1**((p-q-1)::NNI))::UP(x,R)
         List2:L UP(x,R):=append(List2:L UP(x,R),[Pr1]:L UP(x,R))
       List2


-- Computation of the sign (+1,0,-1) of an element in an ordered integral
-- domain

--     sign(r:R):R ==
--       r =$R 0 => 0
--       r >$R 0 => 1
--       -1


-- Computation of the delta function:

     delta(int1:NNI):R ==
       (-1)**((int1*(int1+1) exquo 2)::NNI)


-- Computation of the Sturm-Habicht sequence of two polynomials P and Q
-- in R[x] where R is an ordered integral domaine

     polsth1(p1,p:NNI,p2,q:NNI,c1:R):L UP(x,R) ==
       sc1:R:=(sign(c1))::R
       Pr1:UP(x,R):=pseudoRemainder(differentiate(p1)*p2,p1)
       Pr2:UP(x,R):=(Pr1 exquo c1**(q::NNI))::UP(x,R)
       c2:R:=leadingCoefficient(Pr2)
       r:NNI:=degree(Pr2)
       Pr3:UP(x,R):=monomial(sc1**((p-r-1)::NNI),0)*p1
       Pr4:UP(x,R):=monomial(sc1**((p-r-1)::NNI),0)*Pr2
       Listf:L UP(x,R):=[Pr3,Pr4]
       if r < p-1 then
         Pr5:UP(x,R):=monomial(delta((p-r-1)::NNI)*c2**((p-r-1)::NNI),0)*Pr2
         for j in ((r+1)::INT)..((p-2)::INT) repeat
           Listf:L UP(x,R):=append(Listf:L UP(x,R),[0]:L UP(x,R))
         Listf:L UP(x,R):=append(Listf:L UP(x,R),[Pr5]:L UP(x,R))
       if Pr1=0 then List1:L UP(x,R):=Listf
                else List1:L UP(x,R):=subresultantSequence(p1,Pr2)
       List2:L UP(x,R):=[]
       for j in 0..((r-1)::INT) repeat
         Pr6:UP(x,R):=monomial(delta((p-j-1)::NNI),0)*List1.((p-j+1)::NNI)
         List2:L UP(x,R):=append([Pr6]:L UP(x,R),List2:L UP(x,R))
       append(Listf:L UP(x,R),List2:L UP(x,R))

     polsth2(p1,p:NNI,p2,q:NNI,c1:R):L UP(x,R) ==
       sc1:R:=(sign(c1))::R
       Pr1:UP(x,R):=monomial(sc1,0)*p1
       Pr2:UP(x,R):=differentiate(p1)*p2
       Pr3:UP(x,R):=monomial(sc1,0)*Pr2
       Listf:L UP(x,R):=[Pr1,Pr3]
       List1:L UP(x,R):=subresultantSequence(p1,Pr2)
       List2:L UP(x,R):=[]
       for j in 0..((p-2)::INT) repeat
         Pr4:UP(x,R):=monomial(delta((p-j-1)::NNI),0)*List1.((p-j+1)::NNI)
         Pr5:UP(x,R):=(Pr4 exquo c1)::UP(x,R)
         List2:L UP(x,R):=append([Pr5]:L UP(x,R),List2:L UP(x,R))
       append(Listf:L UP(x,R),List2:L UP(x,R))

     polsth3(p1,p:NNI,p2,q:NNI,c1:R):L UP(x,R) ==
       sc1:R:=(sign(c1))::R
       q1:NNI:=(q-1)::NNI
       v:NNI:=(p+q1)::NNI
       Pr1:UP(x,R):=monomial(delta(q1::NNI)*sc1**((q+1)::NNI),0)*p1
       Listf:L UP(x,R):=[Pr1]
       List1:L UP(x,R):=subresultantSequence(differentiate(p1)*p2,p1)
       List2:L UP(x,R):=[]
       for j in 0..((p-1)::NNI) repeat
         Pr2:UP(x,R):=monomial(delta((v-j)::NNI),0)*List1.((v-j+1)::NNI)
         Pr3:UP(x,R):=(Pr2 exquo c1)::UP(x,R)
         List2:L UP(x,R):=append([Pr3]:L UP(x,R),List2:L UP(x,R))
       append(Listf:L UP(x,R),List2:L UP(x,R))

     SturmHabichtSequence(p1,p2):L UP(x,R) ==
       p:NNI:=degree(p1)
       q:NNI:=degree(p2)
       c1:R:=leadingCoefficient(p1)
       c1 = 1 or q = 1 => polsth1(p1,p,p2,q,c1)
       q = 0 => polsth2(p1,p,p2,q,c1)
       polsth3(p1,p,p2,q,c1)


-- Computation of the Sturm-Habicht principal coefficients of two
-- polynomials P and Q in R[x] where R is an ordered integral domain

     SturmHabichtCoefficients(p1,p2):L R ==
       List1:L UP(x,R):=SturmHabichtSequence(p1,p2)
--       List2:L R:=[]
       qp:NNI:=#(List1)::NNI
       [coefficient(p,(qp-j)::NNI) for p in List1 for j in 1..qp]
--       for j in 1..qp repeat
--         Ply:R:=coefficient(List1.j,(qp-j)::NNI)
--         List2:L R:=append(List2,[Ply])
--       List2


-- Computation of the number of sign variations of a list of non zero
-- elements in an ordered integral domain

     variation(Lsig:L R):INT ==
       #Lsig = 1 => 0
       elt1:R:=first Lsig
       elt2:R:=Lsig.2
       sig1:R:=(sign(elt1*elt2))::R
       List1:L R:=rest Lsig
       sig1 = 1 => variation List1
       1+variation List1


-- Computation of the number of sign permanences of a list of non zero
-- elements in an ordered integral domain

     permanence(Lsig:L R):INT ==
       #Lsig = 1 => 0
       elt1:R:=first Lsig
       elt2:R:=Lsig.2
       sig1:R:=(sign(elt1*elt2))::R
       List1:L R:=rest Lsig
       sig1 = -1 => permanence List1
       1+permanence List1


-- Computation of the functional W which works over a list of elements
-- in an ordered integral domain, with non zero first element

     qzeros(Lsig:L R):L R ==
       while last Lsig = 0 repeat
         Lsig:L R:=reverse rest reverse Lsig
       Lsig

     epsil(int1:NNI,elt1:R,elt2:R):INT ==
       int1 = 0 => 0
       odd? int1 => 0
       ct1:INT:=if positive? elt1 then 1 else -1
       ct2:INT:=if positive? elt2 then 1 else -1
       ct3:NNI:=(int1 exquo 2)::NNI
       ct4:INT:=(ct1*ct2)::INT
       ((-1)**(ct3::NNI))*ct4

     numbnce(Lsig:L R):NNI ==
       null Lsig => 0
       eltp:R:=Lsig.1
       eltp = 0 => 0
       1 + numbnce(rest Lsig)

     numbce(Lsig:L R):NNI ==
       null Lsig => 0
       eltp:R:=Lsig.1
       not(eltp = 0) => 0
       1 + numbce(rest Lsig)

     wfunctaux(Lsig:L R):INT ==
       null Lsig => 0
       List2:L R:=[]
       List1:L R:=Lsig:L R
       cont1:NNI:=numbnce(List1:L R)
       for j in 1..cont1 repeat
         List2:L R:=append(List2:L R,[first List1]:L R)
         List1:L R:=rest List1
       ind2:INT:=0
       cont2:NNI:=numbce(List1:L R)
       for j in 1..cont2 repeat
         List1:L R:=rest List1
         ind2:INT:=epsil(cont2:NNI,last List2,first List1)
       ind3:INT:=permanence(List2:L R)-variation(List2:L R)
       ind4:INT:=ind2+ind3
       ind4+wfunctaux(List1:L R)

     wfunct(Lsig:L R):INT ==
       List1:L R:=qzeros(Lsig:L R)
       wfunctaux(List1:L R)


-- Computation of the integer number:
--    #[{a in Rc(R)/P(a)=0 Q(a)>0}] - #[{a in Rc(R)/P(a)=0 Q(a)<0}]
-- where:
--    - R is an ordered integral domain,
--    - Rc(R) is the real clousure of R,
--    - P and Q are polynomials in R[x],
--    - by #[A] we note the cardinal of the set A

-- In particular:
--      - SturmHabicht(P,1) is the number of "real" roots of P,
--      - SturmHabicht(P,Q**2) is the number of "real" roots of P making Q neq 0

     SturmHabicht(p1,p2):INT ==
--     print("+" :: Ex)
       p2 = 0 => 0
       degree(p1:UP(x,R)) = 0 => 0
       List1:L UP(x,R):=SturmHabichtSequence(p1,p2)
       qp:NNI:=#(List1)::NNI
       wfunct [coefficient(p,(qp-j)::NNI) for p in List1 for j in 1..qp]

     countRealRoots(p1):INT == SturmHabicht(p1,1)

     if R has GcdDomain then
        SturmHabichtMultiple(p1,p2):INT ==
   --     print("+" :: Ex)
          p2 = 0 => 0
          degree(p1:UP(x,R)) = 0 => 0
          SH:L UP(x,R):=SturmHabichtSequence(p1,p2)
          qp:NNI:=#(SH)::NNI
          ans:= wfunct [coefficient(p,(qp-j)::NNI) for p in SH for j in 1..qp]
          SH:=reverse SH
          while first SH = 0 repeat SH:=rest SH
          degree first SH = 0 => ans
          -- OK: it probably wasn't square free, so this item is probably the 
          -- gcd of p1 and p1'
          -- unless p1 and p2 have a factor in common (naughty!)
          differentiate(p1) exquo first SH case UP(x,R) =>
             -- it was the gcd of p1 and p1'
             ans+SturmHabichtMultiple(first SH,p2)
          sqfr:=factorList squareFree p1
          #sqfr = 1 and sqfr.first.xpnt=1 => ans
          reduce("+",[f.xpnt*SturmHabicht(f.fctr,p2) for f in sqfr])

        countRealRootsMultiple(p1):INT == SturmHabichtMultiple(p1,1)

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