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

)abbrev package NORMPK NormalizationPackage
++ Author: Marc Moreno Maza
++ Date Created: 09/23/1998
++ Date Last Updated: 12/16/1998
++ References:
++  [1] D. LAZARD "A new method for solving algebraic systems of 
++      positive dimension" Discr. App. Math. 33:147-160,1991
++  [2] M. MORENO MAZA and R. RIOBOO "Computations of gcd over
++      algebraic towers of simple extensions" In proceedings of AAECC11
++      Paris, 1995.
++  [3] M. MORENO MAZA "Calculs de pgcd au-dessus des tours
++      d'extensions simples et resolution des systemes d'equations
++      algebriques" These, Universite P.etM. Curie, Paris, 1997.
++ Description: 
++ A package for computing normalized assocites of univariate polynomials
++ with coefficients in a tower of simple extensions of a field.

NormalizationPackage(R,E,V,P,TS) : SIG == CODE where
  R : GcdDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R,E,V)
  TS : RegularTriangularSetCategory(R,E,V,P)

  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  S ==> String
  K ==> Fraction R
  LP ==> List P
  PWT ==> Record(val : P, tower : TS)

  BWT ==> Record(val : Boolean, tower : TS)
  LpWT ==> Record(val : (List P), tower : TS)
  Split ==> List TS
  polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P)
  regsetgcdpack ==> SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,TS)

  SIG ==> with

     recip : (P, TS) -> Record(num:P,den:P)
       ++ \axiom{recip(p,ts)} returns the inverse of \axiom{p} w.r.t \spad{ts}
       ++ assuming that \axiom{p} is invertible w.r.t \spad{ts}.

     normalizedAssociate : (P, TS) -> P
       ++ \axiom{normalizedAssociate(p,ts)} returns a normalized polynomial 
       ++ \axiom{n} w.r.t. \spad{ts} such that \axiom{n} and \axiom{p} are
       ++ associates w.r.t \spad{ts} and assuming that \axiom{p} is invertible 
       ++ w.r.t \spad{ts}.

     normalize : (P, TS) -> List PWT
       ++ \axiom{normalize(p,ts)} normalizes \axiom{p} w.r.t \spad{ts}.

     outputArgs : (S, S, P, TS) -> Void
       ++ \axiom{outputArgs(s1,s2,p,ts)} 
       ++ is an internal subroutine, exported only for developement.

     normInvertible? : (P, TS) -> List BWT
       ++ \axiom{normInvertible?(p,ts)} 
       ++ is an internal subroutine, exported only for developement.

  CODE ==> add

     if TS has SquareFreeRegularTriangularSetCategory(R,E,V,P)
     then

       normInvertible?(p:P, ts:TS): List BWT ==
         stoseInvertible?_sqfreg(p,ts)$regsetgcdpack

     else

       normInvertible?(p:P, ts:TS): List BWT ==
         stoseInvertible?_reg(p,ts)$regsetgcdpack

     if (R has RetractableTo(Integer)) and (V has ConvertibleTo(Symbol))
     then 

       outputArgs(s1:S, s2: S, p:P,ts:TS): Void ==
         if not empty? s1 then output(s1, p::OutputForm)$OutputPackage
         if not empty? s1 then _
              output(s1,(convert(p)@String)::OutputForm)$OutputPackage
         output(" ")$OutputPackage
         if not empty? s2 then output(s2, ts::OutputForm)$OutputPackage       
         empty? s2 => void()
         output(s2,("[")::OutputForm)$OutputPackage
         lp: List P := members(ts)
         for q in lp repeat
            output((convert(q)@String)::OutputForm)$OutputPackage
         output("]")$OutputPackage
         output(" ")$OutputPackage

     else

       outputArgs(s1:S, s2: S, p:P,ts:TS): Void ==
         if not empty? s1 then output(s1, p::OutputForm)$OutputPackage
         output(" ")$OutputPackage
         if not empty? s2 then output(s2, ts::OutputForm)$OutputPackage       
         output(" ")$OutputPackage

     recip(p:P,ts:TS): Record(num:P, den:P) ==
     -- ASSUME p is invertible w.r.t. ts
     -- ASSUME mvar(p) is algebraic w.r.t. ts
       v := mvar(p)
       ts_v := select(ts,v)::P
       if mdeg(p) < mdeg(ts_v)
         then
           hesrg: Record (gcd : P, coef2 : P)  := _
                    halfExtendedSubResultantGcd2(ts_v,p)$P
           d: P :=  hesrg.gcd; n: P := hesrg.coef2
         else
           hesrg: Record (gcd : P, coef1 : P) := _
                    halfExtendedSubResultantGcd1(p,ts_v)$P
           d: P :=  hesrg.gcd; n: P := hesrg.coef1
       g := gcd(n,d)
       (n, d) := ((n exquo g)::P, (d exquo g)::P)
       remn, remd: Record(rnum:R,polnum:P,den:R)
       remn := remainder(n,ts); remd := remainder(d,ts)
       cn := remn.rnum; pn := remn.polnum; dn := remn.den
       cd := remd.rnum; pd := remd.polnum; dp := remd.den
       k: K := (cn / cd) * (dp / dn)
       pn := removeZero(pn,ts)
       pd := removeZero(pd,ts)
       [numer(k) * pn, denom(k) * pd]$Record(num:P, den:P)

     normalizedAssociate(p:P,ts:TS): P ==
     -- ASSUME p is invertible or zero w.r.t. ts
       empty? ts => p
       zero?(p) => p
       ground?(p) => 1
       zero? initiallyReduce(init(p),ts) =>
         error "in normalizedAssociate$NORMPK: bad #1"
       vp := mvar(p)
       ip: P := p
       mp: P := 1
       tp: P := 0
       while not ground?(ip) repeat
         v := mvar(ip)
         if algebraic?(v,ts)
           then
             if v = vp
               then
                 ts_v := select(ts,v)::P
                 ip := lastSubResultant(ip,ts_v)$P
                 ip := remainder(ip,ts).polnum
                 -- ip := primitivePart stronglyReduce(ip,ts)
                 ip := primitivePart initiallyReduce(ip,ts)
               else
                 qr := recip(ip,ts)
                 ip := qr.den
                 tp := qr.num * tp
                 zero? ip =>
                     outputArgs("p = ", " ts = ",p,ts)
                     error _
                       "in normalizedAssociate$NORMPK: should never happen !"
           else
             tp := tail(ip) * mp + tp
             mp := mainMonomial(ip) * mp
             ip := init(ip)
       r := ip * mp + tp
       r := remainder(r,ts).polnum
       -- primitivePart stronglyReduce(r,ts)
       primitivePart initiallyReduce(r,ts)

     normalize(p: P, ts: TS): List PWT ==
       zero? p => [[p,ts]$PWT]
       ground? p => [[1,ts]$PWT]
       zero? initiallyReduce(init(p),ts) =>
         error "in normalize$NORMPK: init(#1) reduces to 0 w.r.t. #2"
       --output("Entering  normalize")$OutputPackage
       --outputArgs("p = ", " ts = ",p,ts)
       --output("Calling  normInvertible?")$OutputPackage
       lbwt: List BWT := normInvertible?(p,ts)
       --output("Result is: ")$OutputPackage
       --output(lbwt::OutputForm)$OutputPackage
       lpwt: List PWT := []
       for bwt in lbwt repeat
         us := bwt.tower
         q := remainder(p,us).polnum
         q := removeZero(q,us)
         bwt.val =>
           --output("Calling  normalizedAssociate")$OutputPackage
           --outputArgs("q = ", " us = ",q,us)
           lpwt := cons([normalizedAssociate(q,us)@P,us]$PWT, lpwt)
           --output("Leaving  normalizedAssociate")$OutputPackage
         zero? q => lpwt := cons([0$P,us]$PWT, lpwt)
         lpwt := concat(normalize(q,us)@(List PWT),lpwt)
       lpwt