)ab domain ROINT RightOpenInterval


-- RightOpenInterval
++ Author: Renaud Rioboo
++ Date Created: Summer 1991
++ Date Last Updated: Fri Apr 17 1992
++ Related Constructors: SEGMENT
++ Also See:
++ AMS Classifications:
++ Keywords: Intervals,
++ Basic Functions: provides usual intervals.
++ References: 
++ Description: This domain implements intervals in the usual sense

RightOpenInterval(TheField) : PUB == PRIV where

  TheField : Join(OrderedRing,Field)
  E           ==> OutputForm
  N           ==> NonNegativeInteger
  Z           ==> Integer
  O           ==> OutputForm
  B           ==> Boolean
  SEG         ==> Segment TheField

  PUB == Join(SetCategory,
              CoercibleTo(Segment(TheField))) with
      lo :         $  -> TheField
      hi :         $  -> TheField
      coerce :    SEG ->    $     
      empty?:      $  ->    B
      map:  (TheField->TheField,$) -> $
      SEGMENT : (TheField,TheField) -> $

  PRIV == add

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

   lo(int) == int.low

   hi(int) == int.high

   empty?(int) == int.low >= int.high

   int1 = int2 == (int1.low = int2.low) and (int1.high = int2.high)

   coerce(int:$):E == commasep
     [ hconcat("[" :: E , (int.low)::E), 
       hconcat((int.high)::E,"[" ::E ) ]

   coerce(int:$):SEG == (int.low)..(int.high)

   coerce(seg:SEG):$ == [lo(seg),hi(seg)]

   map(foo,int) == [foo(int.low),foo(int.high)]

   a..b == [a,b]

)ab domain SOAE SimpleOrderedAlgebraicExtension

-- SimpleOrderedAlgebraicExtension

++ Author: Renaud Rioboo
++ Date Created: summer 1988
++ Date Last Updated: Fri Apr 17 1992
++ Basic Functions: provides computations in an ordered algabraic extension 
++ Related Constructors: SparseUnivariatePolynomial, RightOpenInterval
++ Also See:
++ AMS Classifications:
++ Keywords: Real Algebraic Numbers
++ References: 
++ Description:
++ This domain implements an Ordered Extension of an ordered field.
++ The base field is assumed to be archimedian
++ Note: 
++ seg is assumed to be such that pol(hi(seg)) ^= 0
++ The defining polynomial is updated this is implementation dependant
 

SimpleOrderedAlgebraicExtension(TheField,seg,pol,extvar) : PUB == PRIV where

  SEG         ==> RightOpenInterval(TheField)
  E           ==> OutputForm
  P           ==> SparseUnivariatePolynomial(TheField)
  TheField : Join(OrderedRing,Field)
  seg : SEG
  pol : P
  extvar : E
  N           ==> NonNegativeInteger
  Z           ==> Integer
  PINT        ==> List Record(k:N,c:TheField)
  SINT        ==> Record(le:TheField,ri:TheField)
  O           ==> OutputForm
  B           ==> Boolean

  PUB == Join(OrderedRing,
              CommutativeRing,
              ConvertibleTo(P),
              Field,
              FullyRetractableTo(TheField),
              Algebra TheField) with
         poly :                 ()         ->  P
         ++ poly returns the current defining polynomial
         generator :            ()         ->  $
         ++ generator returns the canonical root of poly()
         segment :              ()         -> SEG
         ++ segment returns the current isolating interval of generator()
         changePoly :            P         ->  P
         ++ changePoly updates poly() to a new value
         changeSegment :        SEG        -> SEG
         ++ changeSegment updates segment() to a new value
         refine:                ()         -> SEG
         ++ refine splits segment() into two equal parts
         refine:             TheField      -> SEG
         ++ refine splits segment() using argument as splitting point
         chooseFactor:          P          ->  B
         ++ chooseFactor updates poly() using argument that is a factor 
         ++ chooseFactor returns true if factor is the good one
         starOf :               P          ->  P
         ++ starOf(p)(X) == (X+1)^degree(poly())*p(hi(seg())*X+lo(seg()) / X +1)
         varStar :              P          ->  N
         ++  varStar sign variations of argument's coefficients  
         coerce :               P          ->  $
         ++ coerce to change the type

  PRIV == add
     Rep := P

     coerce(p:P):$ == coerce(p @ $) 

     convert(a) == a pretend P

     characteristic() == 0

     (a:$ * b:$):$ == coerce((a pretend P) *$P (b pretend P)) -- @#$%^ compiler

     (n:Z * a:$):$ == n :: TheField * a

     (rn:TheField * a:$):$ == rn :: $ * a

     coerce(a:$):$ == ((a pretend P) rem$P pol)::$

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

     coerce(a:$):E == outputForm(a ,extvar)

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

     recip(a) ==
        qq:Rep := a 
        d := principalIdeal([pol,qq])$P
        if (degree(d.generator) ^= 0) then
           chooseFactor(d.generator) => "failed"
           d := principalIdeal([pol , qq])$P
        d.coef.2 :: $

     generator() == coerce(monomial(1,1)$P)

     chooseFactor(p) == -- p divides pol that chooses the good one
       degree(p) = 0 => false -- always divide but never factor
       ((p1 := pol exquo$Rep p) case "failed") => error "This ShouldNotAppear"
       -- try to speed the linear case to avoid evaluation on bounds
       degree(p) = 1 =>
             c := - coefficient(p,0) /$TheField leadingCoefficient(p)
             if (^(c < lo(seg)) and (c < hi(seg))) --inside bounds
             then
                changePoly(p)
                true
             else
                changePoly(p1 pretend P)
                false
       degree(p1 pretend P) = 1 =>
             c := - coefficient(p1 pretend P,0) /$TheField 
                  leadingCoefficient(p1 pretend P)
             if (^(c < lo(seg)) and (c < hi(seg))) --inside bounds
             then
                changePoly(p1 pretend P)
                false
             else
                changePoly(p)
                true
       sign((p.lo(seg))*$TheField (p.hi(seg))) <= 0 => 
             changePoly(p)
             true
       changePoly(p1 pretend P)
       false

     a + b == (a +$Rep b)

     -a == (-$Rep a)

     0 == 0$Rep

     1 == 1$Rep

     a = b == zero?(a-b)

     zero?(a) ==
       aa:Rep := a
       degree(aa) = 0 => zero?(leadingCoefficient aa )$TheField
       degree(aa) = 1 => -- to speed up the linear case
           c := - coefficient(aa,0) /$TheField leadingCoefficient(aa)
           if zero?(pol.c)$TheField 
           then         -- root of pol
              chooseFactor(monomial(1,1)-c::P)
           else
              false
       chooseFactor(gcd(pol , aa)$Rep pretend P)


     positive?(a) == positive?(sign(a))$Z

     sign(a) ==
       zero? a => 0$Z
       aa:Rep := a
       degree(aa) = 0 => sign(leadingCoefficient(aa))
       degree(aa) = 1 => -- to speed up the linear case 
          c := - coefficient(aa,0) /$TheField leadingCoefficient(aa)
          while ( ^(c < lo(seg)) and (c < hi(seg)) ) repeat -- inside bounds
              refine()
          sign(aa.lo(seg))
       while (varStar(starOf(aa)) ^= 0) repeat
         refine()
       sign(aa.lo(seg))

     negative?(a) == negative?(sign(a))$Z

     a < b == positive?(b-a)

     varStar(p) == -- if 0 no roots in [0,:infty[
       res : N := 0
       (d := degree(p) @ N) = 0 => res
       lsg := sign(coefficient(p,0))
       for i in 1..d repeat
        if ^( (sg := sign(coefficient(p,i)) ) = 0 ) then
          if (sg ^= lsg) then res := res + 1
          lsg := sg
       res

     starOf(p) == -- P* / P(x) in [a,b[ <=> P*(x) = 0 in [0,+:infty[
       (d := degree(p) @ N) = 0 => p
       n1 :P := monomial(hi(seg),1)$P +$P lo(seg)::P
       n2 :P := 1$P +$P  monomial(1,1)$P
       res : P := coefficient(p,d)::P
       nn2:P := 1$P
       for i in (d-1)..0 by -1 repeat
           nn2 := n2 *$P nn2
           res := n1 *$P res +$P coefficient(p,i pretend N)*$P nn2 
       res



     retractIfCan(a:$):Union(TheField,"failed") == 
        aa : Rep := a
        degree(aa) = 0 => leadingCoefficient(aa)
        degree(pol) = 1 => leadingCoefficient(aa rem$P pol)
        "failed"

     retract(a:$):TheField == 
        (res := retractIfCan(a)) case "failed" => error("Not Retractable")
        res :: TheField

     refine() ==
        refine( (lo(seg)+hi(seg)) /$TheField (1$TheField+1$TheField))

     refine(int) ==
        left := lo(seg)
        right := hi(seg)
        ((int < left) or (int > right)) => error("not in bounds")
        sign(pol.int *$TheField pol.right) <= 0  => changeSegment(int..right)
        changeSegment(left..int)

     poly() == pol

     segment() == seg

     changeSegment(seg1) ==         -- horrible but usefull side effects
       r := seg pretend SINT
       r1 := seg1 pretend SINT      -- implementation dependant
       r.le := r1.le
       r.ri := r1.ri
       seg

     changePoly(pol1) ==
       l := pol pretend PINT        -- horrible but usefull side effects
       l1 := pol1 pretend PINT      -- implementation dependant
       l.first := l1.first
       l.rest := l1.rest
       pol


)ab domain RECLOS RealClosure
-- RealClosure

++ Author: Renaud Rioboo
++ Date Created: summer 1988
++ Date Last Updated: Fri Apr 17 1992
++ Basic Functions: provides computations in an ordered real closure
++ Related Constructors: SimpleOrderedAlgebraicExtension
++ Also See:
++ AMS Classifications:
++ Keywords: Real Algebraic Numbers
++ References: 
++ Description:
++ This domain implements a real closure of an ordered field.
++ The base field is assumed to be archimedian
++ Note: 



RealClosure(TheField): PUB == PRIV where

   TheField : Join(OrderedRing, Field)
   E           ==> OutputForm
   Z           ==> Integer
   SE          ==> Symbol
   P           ==> Polynomial(TheField)
   PME         ==> SparseUnivariatePolynomial($)
   SEG         ==> RightOpenInterval($)
   B           ==> Boolean
   ReRep       ==> SimpleOrderedAlgebraicExtension($,ss,pp,ee)

   -- this is only to avoid writing, segment should be ss, polynomial should
   -- be pp, and variable name should be ee



   PUB == Join(OrderedRing,
               CommutativeRing,
               Field,
               ConvertibleTo PME,
               FullyRetractableTo TheField,
               Algebra TheField) with

       isAlgebraic?   : ($,SE)           -> B
       ++ isAlgebraic? returns true if argument 1 depends on argument 2
       mainVarOf      :   $              -> Union(SE,"failed")
       ++ mainVarOf is the main algebraic quantity of argument
       mainAlgOf      :   $              -> PME
       ++ mainAlgOf is defining poynomial of main algebraic quantity
       mainSegOf      :   $              -> SEG
       ++ mainSegOf is isolating interval of main algebraic quantity
       mainPolOf      :   $              -> Union(PME,"failed")
       ++ mainPolOf isvalue of argument in terms of main algebraic quantity
       lessAlgebraic  : ($,SE)           -> $
       ++ lessAlgebraic removes if it can depency of argument 1 in argument 2
       moreAlgebraic  : (SEG,PME,SE,$)   -> $
       ++ moreAlgebraic return a dummy algebraic element
       algebraicOf    : (SEG,PME,SE)     -> $
       ++ algebraicOf is the generator of a simple ordered extension
       lessAlgebraic  : $                -> $
       ++ lessAlgebraic removes if it can depency of argument in mainAlgOf

   PRIV == add
       Rec := Record(seg: SEG, pol: PME, exp :SE, val:PME)
       Rep := Union(TheField,Rec)


       coerce(x:$):$ ==
          x case TheField => x
          ss := x.seg
          pp := x.pol
          ee := x.exp
          [ss,pp,ee,convert(x.val ::ReRep)]$Rec

       positive?(x) == positive?(sign(x))$Z

       sign(x) ==
          x case TheField => sign(x)$TheField
          ss := x.seg
          pp := x.pol
          ee := x.exp
          sign(x.val :: ReRep)$ReRep          

       negative?(x) == negative?(sign(x))$Z

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

       zero?(x) ==
          x case TheField => x =$TheField 0
          ee := x.exp
          ss := x.seg
          pp := x.pol
          zero?(x.val :: ReRep)$ReRep
 
       x = y == zero?(x-y)

       mainSegOf(x) ==
          x case TheField => (x)..(x)
          x.seg

       mainAlgOf(x) ==
          x case TheField => monomial(1,1)$PME -(x::PME)
          x.pol

       mainVarOf(x) ==
          x case TheField => "failed"
          x.exp

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

       coerce(x:$):E ==
          x case TheField => x::TheField :: E
          xx:$ := coerce(x)
          outputForm(xx.val,x.exp::E)$PME



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


       recip(x) ==
         x case TheField =>
           if ((r := recip(x)$TheField) case TheField)
           then r::$
           else "failed"
         ss := x.seg
         pp := x.pol
         ee := x.exp
         alg := x.val :: ReRep
         if ((r := recip(alg)$ReRep) case "failed")
         then "failed"
         else lessAlgebraic([ss,pp,ee,convert(r::ReRep)]$Rec ::$)

       (n:Z * x:$):$ == 
          (n :: TheField) * x

       (rn:TheField * x:$):$ == 
          rn =$TheField 0 => 0
          rn =$TheField 1 => x
          x case TheField => rn *$TheField x
          -- however wee need to call lessAlgebraic since x.pol can be linear 
          lessAlgebraic([x.seg,x.pol,x.exp,map((rn*#1),x.val)]$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.exp > x.exp => y*moreAlgebraic(y.seg,y.pol,y.exp,x)
          x.exp > y.exp => x*moreAlgebraic(x.seg,x.pol,x.exp,y)
          -- now x.exp = y.exp
          ss := x.seg
          pp := x.pol
          ee := x.exp
          -- we will multiply the polynomials and then transform !
          -- however wee need to call lessAlgebraic since pp can be linear 
          lessAlgebraic([ss,pp,ee, convert((x.val * y.val)::ReRep)]$Rec)

       x + y ==
          (x case TheField) and (y case TheField) => x +$TheField y
          (x case TheField) => 
                lessAlgebraic([y.seg,y.pol,y.exp,x::PME+(y.val)]$Rec)
          -- x is no longer TheField
          (y case TheField) => 
                lessAlgebraic([x.seg,x.pol,x.exp,(x.val)+y::PME]$Rec)
          -- now both are algebraic
          y.exp > x.exp => moreAlgebraic(y.seg,y.pol,y.exp,x)+y
          x.exp > y.exp => x+moreAlgebraic(x.seg,x.pol,x.exp,y)
          -- now x.exp = y.exp
          ss := x.seg
          pp := x.pol
          ee := x.exp
          -- we will add the polynomials and then transform ! 
          -- however wee need to call lessAlgebraic since pp can be linear 
          lessAlgebraic([ss,pp,ee, convert((x.val + y.val)::ReRep)]$Rec)

       -x ==
          x case TheField => -$TheField (x::TheField)
          ss := x.seg
          pp := x.pol
          ee := x.exp
          lessAlgebraic([ss,pp,ee,-(x.val)]$Rec)

       algebraicOf(ss,pp,ee) ==
          [ss, pp, ee, convert(generator()$ReRep)]$Rec


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

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


       isAlgebraic?(x, e) ==
          x case TheField => false
          ee := x.exp
          ee <$SE e => false
          ee >$SE e => _or/[isAlgebraic?(i,e) for i in coefficients(x.val)]
          ss := x.seg
          pp := x.pol
          ee := x.exp
          ^(retractIfCan(x.val::ReRep) case "failed")

       lessAlgebraic(x) ==
          x case TheField => x
          ss := x.seg
          pp := x.pol
          ee := x.exp
          (res:Union($,"failed") := --@#$% compiler
              retractIfCan(x.val::ReRep)) case "failed" => x
          res :: $

       lessAlgebraic(x,e) ==
          x case TheField  => error("Is not Algebraic")
          ss := x.seg
          pp := x.pol
          ee := x.exp
          ee <$SE e => error("Is Not Algebraic")
          ee >$SE e => [ss,pp,ee,map(lessAlgebraic(#1,e), x.val)]$Rec
          (res:Union($,"failed") := --@#$% compiler
              retractIfCan(x.val::ReRep)) case "failed" => 
                    error("Can't Unmake Algebraicity")
          res :: $

       moreAlgebraic(s, p, e, x) ==
          x case TheField => [s, p, e, monomial(x,0)$PME]$Rec
          x.exp <$SE e => [s, p, e, monomial(x,0)$PME]$Rec
          error("Can't make Algebraicity")



       0 == (0$TheField) :: $

       1 == (1$TheField) :: $

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