axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / ROIRC.spad

)abbrev domain ROIRC RightOpenIntervalRootCharacterization
++ Author: Renaud Rioboo
++ Date Created: summer 1992
++ Date Last Updated: January 2004
++ Description:
++ \axiomType{RightOpenIntervalRootCharacterization} provides work with
++ interval root coding.

RightOpenIntervalRootCharacterization(TheField,ThePolDom) : SIG == CODE 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)

  SIG ==> 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
      ++ \axiom{size(rootChar)} is 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

  CODE ==> 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
   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)
     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(z +-> (d-z)::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 o > 0
    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 f*term < 0 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
      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
      q := invert(q)
      hasVars(addOne(q))

   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 (res.high * res.low < 0 ) 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 degree(delta) > 0
     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 (s1*s2 < 0)
        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 s < 0 
        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)