-----------------------------------------------------------------------------
-- 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.
-----------------------------------------------------------------------------
)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)
    
)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

)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))

)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) == 
    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
)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 :: $

