path: root/src/algebra/reclos.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra reclos.spad}
\author{Renaud Rioboo}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
This file describes the Real Closure 1.0 package which consists of different
packages, categoris and domains :

- the package RealPolynomialUtilitiesPackage whichs receives a field and a
univariate polynomial domain with coefficients in the field. It computes some
simple functions such as Strum and Sylvester sequences.

- The category RealRootCharacterizationCategory provides abstarct
functionalities to work with "real roots" of univariate polynomials. These
resemble variables with some functionalities needed to compute important
operations.

- RealClosedField is a category with provides comon operations available over
real closed fiels. These include finding all the roots of univariate
polynomial, taking square roots, ...

- The domain RightOpenIntervalRootCharacterization is the main code that
provides the functionalities of RealRootCharacterizationCategory for the case
of archimedean fileds. Abstract roots are encoded with a left closed right
open interval containing the root together with a defining polynomial for the
root.

- The RealClosure domain is the end-user code, it provides usual arithmetics
with real algebraic numbers, along with the functionalities of a real closed
field. It also provides functions to approximate a real algebraic number by an
element of the base field. This approximation may either be absolute
(approximate) or relative (realtivApprox).


CAVEEATS

Since real algebraic expressions are stored as depending on "real roots" which
are managed like variables, there is an ordering on these. This ordering is
dynamical in the sense that any new algebraic takes precedence over older
ones. In particular every cretaion function raises a new "real root". This has
the effect that when you type something like sqrt(2) + sqrt(2) you have two
new variables which happen to be equal. To avoid this name the expression such
as in s2 := sqrt(2) ; s2 + s2

Also note that computing times depend strongly on the ordering you implicitly
provide. Please provide algebraics in the order which most natural to you.

LIMITATIONS

The file reclos.input show some basic use of the package.  This packages uses
algorithms which are published in [1] and [2] which are based on field
arithmetics, inparticular for polynomial gcd related algorithms. This can be
quite slow for high degree polynomials and subresultants methods usually work
best. Betas versions of the package try to use these techniques in a better
way and work significantly faster. These are mostly based on unpublished
algorithms and cannot be distributed. Please contact the author if you have a
particular problem to solve or want to use these versions.

Be aware that approximations behave as post-processing and that all
computations are done excatly. They can thus be quite time consuming when
depending on several "real roots".
\section{package POLUTIL RealPolynomialUtilitiesPackage}
<<package POLUTIL RealPolynomialUtilitiesPackage>>=
)abbrev package POLUTIL RealPolynomialUtilitiesPackage
++ Author: Renaud Rioboo
++ Date Created: summer 1992
++ Basic Functions: provides polynomial utilities
++ Related Constructors: RealClosure, 
++ Date Last Updated: July 2004
++ Also See: 
++ AMS Classifications:
++ Keywords: Sturm sequences
++ References:  
++ Description:
++ \axiomType{RealPolynomialUtilitiesPackage} provides common functions used
++ by interval coding.
RealPolynomialUtilitiesPackage(TheField,ThePols) : PUB == PRIV where

    TheField : Field
    ThePols : UnivariatePolynomialCategory(TheField)

    Z ==> Integer
    N ==> NonNegativeInteger
    P ==> ThePols

    PUB == with

       sylvesterSequence : (ThePols,ThePols) -> List ThePols
         ++ \axiom{sylvesterSequence(p,q)} is the negated remainder sequence
         ++ of p and q divided by the last computed term
       sturmSequence : ThePols -> List ThePols
         ++ \axiom{sturmSequence(p) = sylvesterSequence(p,p')}
       if TheField has OrderedRing then
         boundOfCauchy : ThePols -> TheField
           ++ \axiom{boundOfCauchy(p)} bounds the roots of p
         sturmVariationsOf : List TheField -> N
           ++ \axiom{sturmVariationsOf(l)} is the number of sign variations 
           ++ in the list of numbers l,
           ++ note that the first term counts as a sign
         lazyVariations : (List(TheField), Z, Z) -> N
           ++ \axiom{lazyVariations(l,s1,sn)} is the number of sign variations 
           ++ in the list of non null numbers [s1::l]@sn,


    PRIV == add

     sturmSequence(p) ==
       sylvesterSequence(p,differentiate(p))

     sylvesterSequence(p1,p2) ==
       res : List(ThePols) := [p1]
       while (p2 ~= 0) repeat
         res := cons(p2 , res)
         (p1 , p2) := (p2 , -(p1 rem p2))
       if positive? degree(p1)
       then
         p1 := unitCanonical(p1)
         res := [ term quo p1 for term in res ]
       reverse! res

     if TheField has OrderedRing
     then

       boundOfCauchy(p) ==
         c :TheField := inv(leadingCoefficient(p))
         l := [ c*term for term in rest(coefficients(p))]
         null(l) => 1
         1 + ("max" / [ abs(t) for t in l ])

--       sturmVariationsOf(l) == 
--         res : N := 0
--         lsg := sign(first(l))
--         for term in l repeat
--           if ^( (sg := sign(term) ) = 0 ) then
--             if (sg ~= lsg) then res := res + 1
--             lsg := sg
--         res

       sturmVariationsOf(l) == 
         null(l) => error "POLUTIL: sturmVariationsOf: empty list !"
         l1 := first(l)
         -- first 0 counts as a sign
         ll : List(TheField) := []
         for term in rest(l) repeat
           -- zeros don't count
           if not(zero?(term)) then ll := cons(term,ll)
         -- if l1 is not zero then ll = reverse(l)
         null(ll) => error "POLUTIL: sturmVariationsOf: Bad sequence"
         ln := first(ll)
         ll := reverse(rest(ll))
         -- if l1 is not zero then first(l) = first(ll)
         -- if l1 is zero then first zero should count as a sign
         zero?(l1) => 1 + lazyVariations(rest(ll),sign(first(ll)),sign(ln))
         lazyVariations(ll, sign(l1), sign(ln))

       lazyVariations(l,sl,sh) ==
         zero?(sl) or zero?(sh) => error "POLUTIL: lazyVariations: zero sign!"
         null(l) =>
           if sl = sh then 0 else 1
         null(rest(l)) => 
           if zero?(first(l))
           then error "POLUTIL: lazyVariations: zero sign!"
           else
             if sl = sh 
             then 
               if (sl = sign(first(l)))
               then 0
               else 2
             -- in this case we save one test
             else 1
         s := sign(l.2)
         lazyVariations([first(l)],sl,s) + 
           lazyVariations(rest(rest(l)),s,sh)
    
@
\section{category RRCC RealRootCharacterizationCategory}
<<category RRCC RealRootCharacterizationCategory>>=
)abbrev category RRCC RealRootCharacterizationCategory
++ Author: Renaud Rioboo
++ Date Created: summer 1992
++ Date Last Updated: January 2004
++ Basic Functions: provides operations with generic real roots of 
++                  polynomials 
++ Related Constructors: RealClosure, RightOpenIntervalRootCharacterization
++ Also See: 
++ AMS Classifications:
++ Keywords: Real Algebraic Numbers
++ References: 
++ Description:
++ \axiomType{RealRootCharacterizationCategory} provides common acces
++ functions for all real root codings.
RealRootCharacterizationCategory(TheField, ThePols ) : Category == PUB where

   TheField : Join(OrderedRing, Field)
   ThePols : UnivariatePolynomialCategory(TheField)

   Z ==> Integer
   N ==> PositiveInteger

   PUB ==>
     SetCategory with

        sign:                ( ThePols, $ )   ->            Z
              ++ \axiom{sign(pol,aRoot)} gives the sign of \axiom{pol}
              ++ interpreted as \axiom{aRoot}
        zero? :              ( ThePols, $ )   ->         Boolean
              ++ \axiom{zero?(pol,aRoot)} answers if \axiom{pol}
              ++ interpreted as \axiom{aRoot} is \axiom{0}
        negative?:           ( ThePols, $ )   ->         Boolean
              ++ \axiom{negative?(pol,aRoot)} answers if \axiom{pol}
              ++ interpreted as \axiom{aRoot} is negative
        positive?:           ( ThePols, $ )   ->         Boolean
              ++ \axiom{positive?(pol,aRoot)} answers if \axiom{pol}
              ++ interpreted as \axiom{aRoot} is positive
        recip:               ( ThePols, $ )   ->   Union(ThePols,"failed") 
              ++ \axiom{recip(pol,aRoot)} tries to inverse \axiom{pol}
              ++ interpreted as \axiom{aRoot}
        definingPolynomial:         $         ->         ThePols
              ++ \axiom{definingPolynomial(aRoot)} gives a polynomial
              ++ such that \axiom{definingPolynomial(aRoot).aRoot = 0} 
        allRootsOf:              ThePols      ->          List $
              ++ \axiom{allRootsOf(pol)} creates all the roots of \axiom{pol} 
              ++ in the Real Closure, assumed in order.
        rootOf:              ( ThePols, N )   ->      Union($,"failed")
              ++ \axiom{rootOf(pol,n)} gives the nth root for the order of the
              ++ Real Closure
        approximate :  (ThePols,$,TheField)   ->    TheField
              ++ \axiom{approximate(term,root,prec)} gives an approximation 
              ++ of \axiom{term} over \axiom{root} with precision \axiom{prec}

        relativeApprox :  (ThePols,$,TheField)   ->    TheField
              ++ \axiom{approximate(term,root,prec)} gives an approximation 
              ++ of \axiom{term} over \axiom{root} with precision \axiom{prec}

      add

        zero?(toTest, rootChar) == 
          sign(toTest, rootChar) = 0
                
        negative?(toTest, rootChar) == 
          negative? sign(toTest, rootChar)
        
        positive?(toTest, rootChar) == 
          positive? sign(toTest, rootChar)

        rootOf(pol,n) ==
          liste:List($):= allRootsOf(pol)
          # liste > n => "failed"
          liste.n

        recip(toInv,rootChar) ==
          degree(toInv) = 0 => 
            res := recip(leadingCoefficient(toInv))
            if (res case "failed") then "failed" else (res::TheField::ThePols)
          defPol := definingPolynomial(rootChar)
          d := principalIdeal([defPol,toInv])
          zero?(d.generator,rootChar) => "failed"
          if (degree(d.generator) ~= 0 )
          then
            defPol := (defPol exquo (d.generator))::ThePols
            d := principalIdeal([defPol,toInv])
          d.coef.2

@
\section{category RCFIELD RealClosedField}
<<category RCFIELD RealClosedField>>=
)abbrev category RCFIELD RealClosedField
++ Author: Renaud Rioboo
++ Date Created: may 1993
++ Date Last Updated: January 2004
++ Basic Functions: provides computations with generic real roots of 
++                  polynomials 
++ Related Constructors: SimpleOrderedAlgebraicExtension, RealClosure
++ Also See: 
++ AMS Classifications:
++ Keywords: Real Algebraic Numbers
++ References: 
++ Description:
++ \axiomType{RealClosedField} provides common acces
++ functions for all real closed fields.
RealClosedField : Category == PUB where

    E ==> OutputForm
    SUP ==> SparseUnivariatePolynomial
    OFIELD ==> Join(OrderedRing,Field)
    PME ==> SUP($)
    N ==> NonNegativeInteger
    PI ==> PositiveInteger
    RN ==> Fraction(Integer)
    Z  ==> Integer
    POLY ==> Polynomial
    PACK ==> SparseUnivariatePolynomialFunctions2

    PUB == Join(CharacteristicZero,
                OrderedRing,
                CommutativeRing,
                Field,
                FullyRetractableTo(Fraction(Integer)),
                Algebra Integer,
                Algebra(Fraction(Integer)),
                RadicalCategory) with

        mainForm :   $ -> Union(E,"failed")
             ++ \axiom{mainForm(x)} is the main algebraic quantity name of 
             ++ \axiom{x}

        mainDefiningPolynomial :   $ -> Union(PME,"failed")
             ++ \axiom{mainDefiningPolynomial(x)} is the defining 
             ++ polynomial for the main algebraic quantity of \axiom{x}

        mainValue :   $ -> Union(PME,"failed")
             ++ \axiom{mainValue(x)} is the expression of \axiom{x} in terms
             ++ of \axiom{SparseUnivariatePolynomial($)} 

        rootOf:          (PME,PI,E)           -> Union($,"failed")
             ++ \axiom{rootOf(pol,n,name)} creates the nth root for the order
             ++ of \axiom{pol} and names it \axiom{name}

        rootOf:          (PME,PI)             -> Union($,"failed")
             ++ \axiom{rootOf(pol,n)} creates the nth root for the order
             ++ of \axiom{pol} and gives it unique name

        allRootsOf:       PME                ->  List $
             ++ \axiom{allRootsOf(pol)} creates all the roots
             ++ of \axiom{pol} naming each uniquely

        allRootsOf:       (SUP(RN))          ->  List $
             ++ \axiom{allRootsOf(pol)} creates all the roots
             ++ of \axiom{pol} naming each uniquely

        allRootsOf:       (SUP(Z))          ->  List $
             ++ \axiom{allRootsOf(pol)} creates all the roots
             ++ of \axiom{pol} naming each uniquely

        allRootsOf:       (POLY($))         ->  List $
             ++ \axiom{allRootsOf(pol)} creates all the roots
             ++ of \axiom{pol} naming each uniquely

        allRootsOf:       (POLY(RN))        ->  List $
             ++ \axiom{allRootsOf(pol)} creates all the roots
             ++ of \axiom{pol} naming each uniquely

        allRootsOf:       (POLY(Z))         ->  List $
             ++ \axiom{allRootsOf(pol)} creates all the roots
             ++ of \axiom{pol} naming each uniquely

        sqrt:            ($,PI)                ->     $
             ++ \axiom{sqrt(x,n)} is \axiom{x ** (1/n)}

        sqrt:              $                  ->     $
             ++ \axiom{sqrt(x)} is \axiom{x ** (1/2)}

        sqrt:             RN                  ->     $
             ++ \axiom{sqrt(x)} is \axiom{x ** (1/2)}

        sqrt:              Z                  ->     $
             ++ \axiom{sqrt(x)} is \axiom{x ** (1/2)}

        rename! :        ($,E)                ->     $
             ++ \axiom{rename!(x,name)} changes the way \axiom{x} is printed

        rename :         ($,E)                ->     $
             ++ \axiom{rename(x,name)} gives a new number that prints as name

        approximate:       ($,$) -> RN
              ++ \axiom{approximate(n,p)} gives an approximation of \axiom{n}
              ++ that has precision \axiom{p}

      add

        sqrt(a:$):$ == sqrt(a,2)

        sqrt(a:RN):$ == sqrt(a::$,2)

        sqrt(a:Z):$ == sqrt(a::$,2)

        characteristic == 0

        rootOf(pol,n,o) == 
          r := rootOf(pol,n)
          r case "failed" => "failed"
          rename!(r,o)

        rootOf(pol,n) ==
          liste:List($):= allRootsOf(pol)
          # liste > n => "failed"
          liste.n


        sqrt(x,n) ==
          n = 1 => x
          zero?(x) => 0
          one?(x) => 1 
          if odd?(n)
          then
            r := rootOf(monomial(1,n) - (x :: PME), 1)
          else
            r := rootOf(monomial(1,n) - (x :: PME), 2)
          r case "failed" => error "no roots"
          n = 2 => rename(r,root(x::E)$E)
          rename(r,root(x :: E, n :: E)$E)

        (x : $) ** (rn : RN) == sqrt(x**numer(rn),denom(rn)::PI)

        nthRoot(x, n) == 
          zero?(n) => x
          negative?(n) => inv(sqrt(x,(-n) :: PI))
          sqrt(x,n :: PI)

        allRootsOf(p:SUP(RN)) == allRootsOf(map(#1 :: $ ,p)$PACK(RN,$))

        allRootsOf(p:SUP(Z)) == allRootsOf(map(#1 :: $ ,p)$PACK(Z,$))

        allRootsOf(p:POLY($)) == allRootsOf(univariate(p))

        allRootsOf(p:POLY(RN)) == allRootsOf(univariate(p))

        allRootsOf(p:POLY(Z)) == allRootsOf(univariate(p))

@
\section{domain ROIRC RightOpenIntervalRootCharacterization}
\subsection{makeChar performance problem}
The following lines of code, which check for a possible error,
cause major performance problems and were removed by Renaud Rioboo,
the original author. They were originally inserted for debugging.
\begin{verbatim}
    right <= left => error "ROIRC: makeChar: Bad interval"
    (pol.left * pol.right) > 0 => error "ROIRC: makeChar: Bad pol"
\end{verbatim}
<<performance problem>>=
@
<<domain ROIRC RightOpenIntervalRootCharacterization>>=
)abbrev domain ROIRC RightOpenIntervalRootCharacterization 
++ Author: Renaud Rioboo
++ Date Created: summer 1992
++ Date Last Updated: January 2004
++ Basic Functions: provides computations with real roots of olynomials 
++ Related Constructors: RealRootCharacterizationCategory, RealClosure
++ Also See: 
++ AMS Classifications:
++ Keywords: Real Algebraic Numbers
++ References: 
++ Description:
++ \axiomType{RightOpenIntervalRootCharacterization} provides work with
++ interval root coding.
RightOpenIntervalRootCharacterization(TheField,ThePolDom) : PUB == PRIV where

  TheField : Join(OrderedRing,Field)
  ThePolDom : UnivariatePolynomialCategory(TheField)


  Z           ==>  Integer
  P           ==>  ThePolDom
  N           ==>  NonNegativeInteger
  B           ==>  Boolean
  UTIL        ==>  RealPolynomialUtilitiesPackage(TheField,ThePolDom)
  RRCC        ==>  RealRootCharacterizationCategory
  O ==> OutputForm
  TwoPoints ==> Record(low:TheField , high:TheField)

  PUB == RealRootCharacterizationCategory(TheField, ThePolDom) with

      left    :             $            -> TheField
           ++ \axiom{left(rootChar)} is the left bound of the isolating
           ++ interval
      right   :             $            -> TheField
           ++ \axiom{right(rootChar)} is the right bound of the isolating
           ++ interval
      size    :             $            -> TheField
           ++ The size of the isolating interval
      middle  :             $            -> TheField
           ++ \axiom{middle(rootChar)} is the middle of the isolating
           ++ interval
      refine  :             $            ->    $
           ++ \axiom{refine(rootChar)} shrinks isolating interval around 
           ++ \axiom{rootChar}
      mightHaveRoots :     (P,$)         ->    B
           ++ \axiom{mightHaveRoots(p,r)} is false if \axiom{p.r} is not 0
      relativeApprox :     (P,$,TheField) -> TheField
           ++ \axiom{relativeApprox(exp,c,p) = a} is relatively close to exp
           ++ as a polynomial in c ip to precision p

  PRIV == add



-- local functions


   makeChar:             (TheField,TheField,ThePolDom) ->     $
   refine! :                              $            ->     $
   sturmIsolate : (List(P), TheField, TheField,N,N)    -> List TwoPoints
   isolate :                            List(P)        -> List TwoPoints
   rootBound :                             P           ->   TheField
--   varStar :                                P          ->     N
   linearRecip :                       ( P , $)        -> Union(P, "failed")
   linearZero? :                     (TheField,$)      ->     B
   linearSign :                          (P,$)         ->     Z
   sturmNthRoot : (List(P), TheField, TheField,N,N,N)  -> Union(TwoPoints,"failed")
   addOne :                              P             ->      P
   minus :                               P             ->      P
   translate :                    (P,TheField)         ->      P
   dilate :                       (P,TheField)         ->      P
   invert :                              P             ->      P
   evalOne :                             P             ->   TheField
   hasVarsl:                     List(TheField)        ->      B
   hasVars:                              P             ->      B

-- Representation

   Rep:= Record(low:TheField,high:TheField,defPol:ThePolDom)

-- and now the code !


   size(rootCode) ==
     rootCode.high - rootCode.low

   relativeApprox(pval,rootCode,prec) ==
     -- beurk !
     dPol := rootCode.defPol
     degree(dPol) = 1 => 
       c := -coefficient(dPol,0)/leadingCoefficient(dPol)
       pval.c
     pval := pval rem dPol
     degree(pval) = 0 => leadingCoefficient(pval)
     zero?(pval,rootCode)  => 0
     while mightHaveRoots(pval,rootCode) repeat
          rootCode := refine(rootCode)
     dpval := differentiate(pval)
     degree(dpval) = 0 =>
       l := left(rootCode)
       r := right(rootCode)
       a := pval.l
       b := pval.r
       while ( abs(2*(a-b)/(a+b)) > prec ) repeat
         rootCode := refine(rootCode)
         l := left(rootCode)
         r := right(rootCode)
         a := pval.l
         b := pval.r
       (a+b)/(2::TheField)
     zero?(dpval,rootCode) => 
        relativeApprox(pval, 
                       [left(rootCode),
                         right(rootCode),
                           gcd(dpval,rootCode.defPol)]$Rep,
                       prec)
     while mightHaveRoots(dpval,rootCode) repeat
          rootCode := refine(rootCode)
     l := left(rootCode)
     r := right(rootCode)
     a := pval.l
     b := pval.r
     while ( abs(2*(a-b)/(a+b)) > prec ) repeat
       rootCode := refine(rootCode)
       l := left(rootCode)
       r := right(rootCode)
       a := pval.l
       b := pval.r
     (a+b)/(2::TheField)

   approximate(pval,rootCode,prec) ==
     -- glurp
     dPol := rootCode.defPol
     degree(dPol) = 1 => 
       c := -coefficient(dPol,0)/leadingCoefficient(dPol)
       pval.c
     pval := pval rem dPol
     degree(pval) = 0 => leadingCoefficient(pval)
     dpval := differentiate(pval)
     a : TheField
     b : TheField
     degree(dpval) = 0 =>
       l := left(rootCode)
       r := right(rootCode)
       while ( abs((a := pval.l) - (b := pval.r)) > prec ) repeat
         rootCode := refine(rootCode)
         l := left(rootCode)
         r := right(rootCode)
       (a+b)/(2::TheField)
     zero?(dpval,rootCode) => 
        approximate(pval, 
                    [left(rootCode),
                     right(rootCode),
                      gcd(dpval,rootCode.defPol)]$Rep,
                    prec)
     while mightHaveRoots(dpval,rootCode) repeat
          rootCode := refine(rootCode)
     l := left(rootCode)
     r := right(rootCode)
     while ( abs((a := pval.l) - (b := pval.r)) > prec ) repeat
       rootCode := refine(rootCode)
       l := left(rootCode)
       r := right(rootCode)
     (a+b)/(2::TheField)


   addOne(p) == p.(monomial(1,1)+(1::P))

   minus(p) == p.(monomial(-1,1))

   translate(p,a) == p.(monomial(1,1)+(a::P))

   dilate(p,a) == p.(monomial(a,1))

   evalOne(p) == "+" / coefficients(p)

   invert(p) == 
        d := degree(p)
        mapExponents((d-#1)::N, p)

   rootBound(p) ==
     res : TheField := 1
     raw :TheField := 1+boundOfCauchy(p)$UTIL
     while (res < raw) repeat
       res := 2*(res)
     res

   sturmNthRoot(lp,l,r,vl,vr,n) ==
    nv := (vl - vr)::N
    nv < n => "failed"
    ((nv = 1) and (n = 1)) => [l,r]
    int := (l+r)/(2::TheField)
    lt:List(TheField):=[]
    for t in lp repeat
        lt := cons(t.int , lt)
    vi := sturmVariationsOf(reverse! lt)$UTIL
    o :Z := n - vl + vi
    if positive? o
    then 
       sturmNthRoot(lp,int,r,vi,vr,o::N)
    else
       sturmNthRoot(lp,l,int,vl,vi,n)

   sturmIsolate(lp,l,r,vl,vr) ==
    r <= l => error "ROIRC: sturmIsolate: bad bounds"
    n := (vl - vr)::N
    zero?(n) => []
    one?(n) => [[l,r]]
    int := (l+r)/(2::TheField)
    vi := sturmVariationsOf( [t.int for t in lp ] )$UTIL
    append(sturmIsolate(lp,l,int,vl,vi),sturmIsolate(lp,int,r,vi,vr))

   isolate(lp) ==
     b := rootBound(first(lp))
     l1,l2 : List(TheField)
     (l1,l2) := ([] , [])
     for t in reverse(lp) repeat
       if odd?(degree(t))
       then
        (l1,l2):= (cons(-leadingCoefficient(t),l1),
                   cons(leadingCoefficient(t),l2))
       else
        (l1,l2):= (cons(leadingCoefficient(t),l1),
                   cons(leadingCoefficient(t),l2))
     sturmIsolate(lp,
                  -b,
                  b,
                  sturmVariationsOf(l1)$UTIL,
                  sturmVariationsOf(l2)$UTIL)

   rootOf(pol,n) ==
    ls := sturmSequence(pol)$UTIL
    pol := unitCanonical(first(ls)) -- this one is SqFR
    degree(pol) = 0 => "failed"
    numberOfMonomials(pol) = 1 => ([0,1,monomial(1,1)]$Rep)::$
    b := rootBound(pol)
    l1,l2 : List(TheField)
    (l1,l2) := ([] , [])
    for t in reverse(ls) repeat
      if odd?(degree(t))
      then
       (l1,l2):= (cons(leadingCoefficient(t),l1),
                  cons(-leadingCoefficient(t),l2))
      else
       (l1,l2):= (cons(leadingCoefficient(t),l1),
                  cons(leadingCoefficient(t),l2))
    res := sturmNthRoot(ls,
                        -b,
                        b,
                        sturmVariationsOf(l2)$UTIL,
                        sturmVariationsOf(l1)$UTIL,
                        n)
    res case "failed" => "failed"
    makeChar(res.low,res.high,pol)

   allRootsOf(pol) == 
    ls := sturmSequence(unitCanonical pol)$UTIL
    pol := unitCanonical(first(ls)) -- this one is SqFR
    degree(pol) = 0 => []
    numberOfMonomials(pol) = 1 => [[0,1,monomial(1,1)]$Rep]
    [ makeChar(term.low,term.high,pol) for term in isolate(ls) ]


   hasVarsl(l:List(TheField)) ==
    null(l) => false
    f := sign(first(l))
    for term in rest(l) repeat
      if negative?(f*term) then return(true)
    false
    
   hasVars(p:P) ==
    zero?(p) => error "ROIRC: hasVars: null polynonial"
    zero?(coefficient(p,0)) => true
    hasVarsl(coefficients(p))


   mightHaveRoots(p,rootChar) == 
      a := rootChar.low
      q := translate(p,a)
      not(hasVars(q)) => false
--      varStar(q) = 0 => false
      a := (rootChar.high) - a
      q := dilate(q,a)
      sign(coefficient(q,0))*sign(evalOne(q)) <= 0 => true
      q := minus(addOne(q))
      not(hasVars(q)) => false
--      varStar(q) = 0 => false
      q := invert(q)
      hasVars(addOne(q))
--      ^(varStar(addOne(q)) = 0)

   coerce(rootChar:$):O == 
     commaSeparate([ hconcat("[" :: O , (rootChar.low)::O), 
                     hconcat((rootChar.high)::O,"[" ::O ) ])

   c1 = c2 == 
     mM := max(c1.low,c2.low)
     Mm := min(c1.high,c2.high)
     mM >= Mm => false
     rr : ThePolDom := gcd(c1.defPol,c2.defPol)
     degree(rr) = 0 => false
     sign(rr.mM) * sign(rr.Mm) <= 0

   makeChar(left,right,pol) == 
<<performance problem>>
    res :$ := [left,right,leadingMonomial(pol)+reductum(pol)]$Rep -- safe copy
    while zero?(pol.(res.high)) repeat refine!(res)
    while negative?(res.high * res.low) repeat refine!(res)
    zero?(pol.(res.low)) => [res.low,res.high,monomial(1,1)-(res.low)::P]
    res

   definingPolynomial(rootChar) == rootChar.defPol

   linearRecip(toTest,rootChar) ==
      c := - inv(leadingCoefficient(toTest)) * coefficient(toTest,0)
      r := recip(rootChar.defPol.c)
      if (r case "failed")
      then
        if (c - rootChar.low) * (c - rootChar.high) <= 0
        then 
          "failed"
        else
          newPol := (rootChar.defPol exquo toTest)::P
          ((1$ThePolDom - inv(newPol.c)*newPol) exquo toTest)::P
      else
         ((1$ThePolDom - (r::TheField)*rootChar.defPol) exquo toTest)::P

   recip(toTest,rootChar) ==
     degree(toTest) = 0 or degree(rootChar.defPol) <= degree(toTest) =>
       error "IRC: recip: Not reduced"
     degree(rootChar.defPol) = 1 =>
       error "IRC: recip: Linear Defining Polynomial"
     degree(toTest) = 1 =>
       linearRecip(toTest, rootChar)
     d := extendedEuclidean((rootChar.defPol),toTest)
     (degree(d.generator) = 0 ) => 
         d.coef2
     d.generator := unitCanonical(d.generator)
     (d.generator.(rootChar.low) *
      d.generator.(rootChar.high)<= 0) => "failed"
     newPol := (rootChar.defPol exquo (d.generator))::P
     degree(newPol) = 1 =>
       c := - inv(leadingCoefficient(newPol)) * coefficient(newPol,0)
       inv(toTest.c)::P
     degree(toTest) = 1 => 
       c := - coefficient(toTest,0)/ leadingCoefficient(toTest)
       ((1$ThePolDom - inv(newPol.(c))*newPol) exquo toTest)::P
     d := extendedEuclidean(newPol,toTest)
     d.coef2

   linearSign(toTest,rootChar) ==
      c := - inv(leadingCoefficient(toTest)) * coefficient(toTest,0)
      ev := sign(rootChar.defPol.c)
      if zero?(ev)
      then
        if (c - rootChar.low) * (c - rootChar.high) <= 0
        then
          0
        else
          sign(toTest.(rootChar.high))
      else
        if (ev*sign(rootChar.defPol.(rootChar.high)) <= 0 )
        then
          sign(toTest.(rootChar.high))
        else
          sign(toTest.(rootChar.low))

   sign(toTest,rootChar) ==
     degree(toTest) = 0 or degree(rootChar.defPol) <= degree(toTest) =>
       error "IRC: sign: Not reduced"
     degree(rootChar.defPol) = 1 =>
       error "IRC: sign: Linear Defining Polynomial"
     degree(toTest) = 1 =>
      linearSign(toTest, rootChar)
     s := sign(leadingCoefficient(toTest))
     toTest := monomial(1,degree(toTest))+
               inv(leadingCoefficient(toTest))*reductum(toTest)
     delta := gcd(toTest,rootChar.defPol)
     newChar := [rootChar.low,rootChar.high,rootChar.defPol]$Rep
     if positive? degree(delta)
     then
       if sign(delta.(rootChar.low) * delta.(rootChar.high)) <= 0
       then
        return(0)
       else
        newChar.defPol := (newChar.defPol exquo delta) :: P
        toTest := toTest rem (newChar.defPol)
     degree(toTest) = 0 => s * sign(leadingCoefficient(toTest))
     degree(toTest) = 1 => s * linearSign(toTest, newChar)
     while mightHaveRoots(toTest,newChar) repeat
       newChar := refine(newChar)
     s*sign(toTest.(newChar.low))

   linearZero?(c,rootChar) == 
      zero?((rootChar.defPol).c) and 
       (c - rootChar.low) * (c - rootChar.high) <= 0

   zero?(toTest,rootChar) ==
     degree(toTest) = 0 or degree(rootChar.defPol) <= degree(toTest) =>
       error "IRC: zero?: Not reduced"
     degree(rootChar.defPol) = 1 =>
       error "IRC: zero?: Linear Defining Polynomial"
     degree(toTest) = 1 => 
      linearZero?(- inv(leadingCoefficient(toTest)) * coefficient(toTest,0),
                  rootChar)
     toTest := monomial(1,degree(toTest))+
               inv(leadingCoefficient(toTest))*reductum(toTest)
     delta := gcd(toTest,rootChar.defPol)
     degree(delta) = 0 => false
     sign(delta.(rootChar.low) * delta.(rootChar.high)) <= 0


   refine!(rootChar) ==
     -- this is not a safe function, it can work with badly created object
     -- we do not assume (rootChar.defPol).(rootChar.high) <> 0
        int := middle(rootChar)
        s1 := sign((rootChar.defPol).(rootChar.low))
        zero?(s1) =>
          rootChar.high := int
          rootChar.defPol := monomial(1,1) - (rootChar.low)::P
          rootChar
        s2 := sign((rootChar.defPol).int)
        zero?(s2) =>
          rootChar.low := int
          rootChar.defPol := monomial(1,1) - int::P
          rootChar
        if negative?(s1*s2)
        then 
          rootChar.high := int
        else 
          rootChar.low := int
        rootChar

   refine(rootChar) ==
     -- we assume (rootChar.defPol).(rootChar.high) <> 0
        int := middle(rootChar)
        s:= (rootChar.defPol).int * (rootChar.defPol).(rootChar.high)
        zero?(s) => [int,rootChar.high,monomial(1,1)-int::P]
        if negative? s 
        then 
          [int,rootChar.high,rootChar.defPol]
        else 
          [rootChar.low,int,rootChar.defPol]

   left(rootChar) == rootChar.low

   right(rootChar) == rootChar.high

   middle(rootChar) == (rootChar.low + rootChar.high)/(2::TheField)

--   varStar(p) == -- if 0 no roots in [0,:infty[
--     res : N := 0
--     lsg := sign(coefficient(p,0))
--     l := [ sign(i) for i in reverse!(coefficients(p))]
--     for sg in l repeat
--      if (sg ~= lsg) then res := res + 1
--      lsg := sg
--     res
@
\section{domain RECLOS RealClosure}
The domain constructore {\bf RealClosure} by Renaud Rioboo (University
of Paris 6, France) provides the real closure of an ordered field.
The implementation is based on interval arithmetic. Moreover, the
design of this constructor and its related packages allows an easy
use of other codings for real algebraic numbers.
ordered field
<<domain RECLOS RealClosure>>=
)abbrev domain RECLOS RealClosure
++ Author: Renaud Rioboo
++ Date Created: summer 1988
++ Date Last Updated: January 2004
++ Basic Functions: provides computations in an ordered real closure
++ Related Constructors: RightOpenIntervalRootCharacterization
++ Also See:
++ AMS Classifications:
++ Keywords: Real Algebraic Numbers
++ References: 
++ Description:
++ This domain implements the real closure of an ordered field.
++ Note: 
++ The code here is generic i.e. it does not depend of the way the operations
++ are done. The two macros PME and SEG should be passed as functorial
++ arguments to the domain. It does not help much to write a category
++ since non trivial methods cannot be placed there either.
++ 
RealClosure(TheField): PUB == PRIV where

   TheField   : Join(OrderedRing, Field, RealConstant)

--   ThePols    : UnivariatePolynomialCategory($)
--   PME       ==> ThePols
--   TheCharDom : RealRootCharacterizationCategory($, ThePols )
--   SEG       ==> TheCharDom
-- this does not work yet

   E         ==> OutputForm
   Z         ==> Integer
   SE        ==> Symbol
   B         ==> Boolean
   SUP       ==> SparseUnivariatePolynomial($)
   N         ==> PositiveInteger
   RN        ==> Fraction Z
   LF        ==> ListFunctions2($,N)

-- *****************************************************************
-- *****************************************************************
--             PUT YOUR OWN PREFERENCE HERE
-- *****************************************************************
-- *****************************************************************
   PME       ==> SparseUnivariatePolynomial($)
   SEG       ==> RightOpenIntervalRootCharacterization($,PME)
-- *****************************************************************
-- *****************************************************************


   PUB == Join(RealClosedField,
               FullyRetractableTo TheField,
               Algebra TheField) with

       algebraicOf :   (SEG,E) -> $
             ++ \axiom{algebraicOf(char)} is the external number

       mainCharacterization :   $ -> Union(SEG,"failed")
             ++ \axiom{mainCharacterization(x)} is the main algebraic 
             ++ quantity of \axiom{x} (\axiom{SEG})

       relativeApprox :     ($,$) -> RN
             ++ \axiom{relativeApprox(n,p)} gives a relative 
             ++ approximation of \axiom{n} 
             ++ that has precision \axiom{p}

   PRIV == add

-- local functions

       lessAlgebraic  : $ -> $
       newElementIfneeded : (SEG,E) -> $

-- Representation

       Rec := Record(seg: SEG, val:PME, outForm:E, order:N)
       Rep := Union(TheField,Rec)

-- global (mutable) variables

       orderOfCreation : N := 1$N
          -- it is internally used to sort the algebraic levels

       instanceName : Symbol := new()$Symbol
          -- this used to print the results, thus different instanciations
          -- use different names

-- now the code

       relativeApprox(nbe,prec) ==
          nbe case TheField => retract(nbe)
          appr := relativeApprox(nbe.val, nbe.seg, prec)
          -- now appr has the good exact precision but is $
          relativeApprox(appr,prec)


       approximate(nbe,prec) ==
          abs(nbe) < prec => 0
          nbe case TheField => retract(nbe)
          appr := approximate(nbe.val, nbe.seg, prec)
          -- now appr has the good exact precision but is $
          approximate(appr,prec)

       newElementIfneeded(s,o) ==
         p := definingPolynomial(s)
         degree(p) = 1 => 
             - coefficient(p,0) / leadingCoefficient(p)
         res := [s, monomial(1,1), o, orderOfCreation ]$Rec
         orderOfCreation := orderOfCreation + 1
         res :: $

       algebraicOf(s,o) ==
         pol := definingPolynomial(s)
         degree(pol) = 1 => 
           -coefficient(pol,0) / leadingCoefficient(pol) 
         res := [s, monomial(1,1), o, orderOfCreation ]$Rec
         orderOfCreation := orderOfCreation + 1
         res :: $
         
       rename!(x,o) ==
         x.outForm := o
         x

       rename(x,o) ==
         [x.seg, x.val, o, x.order]$Rec

       rootOf(pol,n) ==
        degree(pol) = 0 => "failed"
        degree(pol) = 1 =>
          if n=1
          then
            -coefficient(pol,0) / leadingCoefficient(pol)
          else
            "failed"
        r := rootOf(pol,n)$SEG
        r case "failed" => "failed"
        o := hconcat(instanceName :: E , orderOfCreation :: E)$E
        algebraicOf(r,o)

       allRootsOf(pol:SUP):List($) == 
        degree(pol)=0 => []
        degree(pol)=1 => [-coefficient(pol,0) / leadingCoefficient(pol)]
        liste := allRootsOf(pol)$SEG
        res : List $ := []
        for term in liste repeat
           o := hconcat(instanceName :: E , orderOfCreation :: E)$E
           res := cons(algebraicOf(term,o), res)
        reverse! res

       coerce(x:$):$ ==
          x case TheField => x
          [x.seg,x.val rem$PME definingPolynomial(x.seg),x.outForm,x.order]$Rec

       positive?(x) == 
          x case TheField => positive?(x)$TheField
          positive?(x.val,x.seg)$SEG

       negative?(x) == 
          x case TheField => negative?(x)$TheField
          negative?(x.val,x.seg)$SEG

       abs(x) == sign(x)*x

       sign(x) ==
          x case TheField => sign(x)$TheField
          sign(x.val,x.seg)$SEG

       x < y == positive?(y-x)

       x = y == zero?(x-y)

       mainCharacterization(x) ==
          x case TheField => "failed"
          x.seg

       mainDefiningPolynomial(x) ==
          x case TheField => "failed"
          definingPolynomial x.seg

       mainForm(x) ==
          x case TheField => "failed"
          x.outForm

       mainValue(x) ==
          x case TheField => "failed"
          x.val

       coerce(x:$):E ==
          x case TheField => x::TheField :: E
          xx:$ := coerce(x)
          outputForm(univariate(xx.val),x.outForm)$SUP


       inv(x) ==
          (res:= recip x) case "failed" => error "Division by 0"
          res :: $

       recip(x) ==
         x case TheField =>
           if ((r := recip(x)$TheField) case TheField)
           then r::$
           else "failed"
         if ((r := recip(x.val,x.seg)$SEG) case "failed")
         then "failed"
         else lessAlgebraic([x.seg,r::PME,x.outForm,x.order]$Rec) 

       (n:Z * x:$):$ == 
          x case TheField => n *$TheField x
          zero?(n) => 0
          one?(n) => x
          [x.seg,map(n * #1, x.val),x.outForm,x.order]$Rec

       (rn:TheField * x:$):$ == 
          x case TheField => rn *$TheField x
          zero?(rn) => 0
          one?(rn) => x
          [x.seg,map(rn * #1, x.val),x.outForm,x.order]$Rec

       (x:$ * y:$):$ ==
          (x case TheField) and (y case TheField) => x *$TheField y
          (x case TheField) => x::TheField * y
              -- x is no longer TheField
          (y case TheField) => y::TheField * x
              -- now both are algebraic
          y.order > x.order => 
            [y.seg,map(x * #1 , y.val),y.outForm,y.order]$Rec
          x.order > y.order => 
            [x.seg,map( #1 * y , x.val),x.outForm,x.order]$Rec
              -- now x.exp = y.exp
              -- we will multiply the polynomials and then reduce
              -- however wee need to call lessAlgebraic  
          lessAlgebraic([x.seg,
                         (x.val * y.val) rem definingPolynomial(x.seg),
                         x.outForm,
                         x.order]$Rec)

       nonNull(r:Rec):$ ==
         degree(r.val)=0 => leadingCoefficient(r.val)
         numberOfMonomials(r.val) = 1 => r
         zero?(r.val,r.seg)$SEG => 0
         r

--       zero?(x) ==
--          x case TheField => zero?(x)$TheField
--          zero?(x.val,x.seg)$SEG
 
       zero?(x) ==
          x case TheField => zero?(x)$TheField
          false
 
       x + y ==
          (x case TheField) and (y case TheField) => x +$TheField y
          (x case TheField) => 
             if zero?(x)
             then 
               y
             else 
               nonNull([y.seg,x::PME+(y.val),y.outForm,y.order]$Rec)
             -- x is no longer TheField
          (y case TheField) => 
             if zero?(y)
             then 
               x
             else 
               nonNull([x.seg,(x.val)+y::PME,x.outForm,x.order]$Rec)
             -- now both are algebraic
          y.order > x.order => 
               nonNull([y.seg,x::PME+y.val,y.outForm,y.order]$Rec)
          x.order > y.order => 
               nonNull([x.seg,(x.val)+y::PME,x.outForm,x.order]$Rec)
              -- now x.exp = y.exp 
              -- we simply add polynomials (since degree cannot increase)
              -- however wee need to call lessAlgebraic  
          nonNull([x.seg,x.val + y.val,x.outForm,x.order])


       -x ==
          x case TheField => -$TheField (x::TheField)
          [x.seg,-$PME x.val,x.outForm,x.order]$Rec


       retractIfCan(x:$):Union(TheField,"failed") ==
          x case TheField => x
          o := x.order
          res := lessAlgebraic x
          res case TheField => res
          o = res.order => "failed"
          retractIfCan res

       retract(x:$):TheField ==
          x case TheField => x
          o := x.order
          res := lessAlgebraic x
          res case TheField => res
          o = res.order => error "Can't retract"
          retract res


       lessAlgebraic(x) ==
          x case TheField => x
          degree(x.val) = 0 => leadingCoefficient(x.val)
          def := definingPolynomial(x.seg)
          degree(def) = 1 => 
            x.val.(- coefficient(def,0) / leadingCoefficient(def))
          x

       0 == (0$TheField) :: $

       1 == (1$TheField) :: $

       coerce(rn:TheField):$ == rn :: $

@
\section{License}
<<license>>=
-----------------------------------------------------------------------------
-- This software was written by Renaud Rioboo (Computer Algebra group of
-- Laboratoire d'Informatique de Paris 6) and is the property of university
-- Paris 6.
-----------------------------------------------------------------------------
@
<<*>>=
<<license>>
<<package POLUTIL RealPolynomialUtilitiesPackage>>
<<category RRCC RealRootCharacterizationCategory>>
<<category RCFIELD RealClosedField>>
<<domain ROIRC RightOpenIntervalRootCharacterization>>
<<domain RECLOS RealClosure>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} R. Rioboo,
{\sl Real Algebraic Closure of an ordered Field : Implementation in Axiom.},
In proceedings of the ISSAC'92 Conference, Berkeley 1992 pp. 206-215.
\bibitem{2} Z. Ligatsikas, R. Rioboo, M. F. Roy 
{\sl Generic computation of the real closure of an ordered field.},
In Mathematics and Computers in Simulation Volume 42, Issue 4-6,
November 1996.
\end{thebibliography}
\end{document}