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

)abbrev domain NSMP NewSparseMultivariatePolynomial
++ Author: Marc Moreno Maza
++ Date Created: 22/04/94
++ Date Last Updated: 14/12/1998
++ Description: 
++ A post-facto extension for \axiomType{SMP} in order
++ to speed up operations related to pseudo-division and gcd.
++ This domain is based on the \axiomType{NSUP} constructor which is 
++ itself a post-facto extension of the \axiomType{SUP} constructor.

NewSparseMultivariatePolynomial(R,VarSet) : SIG == CODE where
  R : Ring
  VarSet : OrderedSet

  N ==> NonNegativeInteger
  Z ==> Integer
  SUP ==> NewSparseUnivariatePolynomial
  SMPR ==> SparseMultivariatePolynomial(R, VarSet)
  SUP2 ==> NewSparseUnivariatePolynomialFunctions2($,$)

  SIG ==> Join(RecursivePolynomialCategory(R,IndexedExponents VarSet,_
                  VarSet), CoercibleTo(SMPR),RetractableTo(SMPR))

  CODE ==> SparseMultivariatePolynomial(R, VarSet) add

     D := NewSparseUnivariatePolynomial($)
     VPoly:=  Record(v:VarSet,ts:D)
     Rep:= Union(R,VPoly)

    --local function
     PSimp: (D,VarSet) -> %

     PSimp(up,mv) ==
       if degree(up) = 0 then leadingCoefficient(up) else [mv,up]$VPoly

     coerce (p:$):SMPR == 
       p pretend SMPR

     coerce (p:SMPR):$ == 
       p pretend $

     retractIfCan (p:$) : Union(SMPR,"failed") == 
       (p pretend SMPR)::Union(SMPR,"failed")

     mvar p == 
       p case R => error" Error in mvar from NSMP : #1 has no variables."
       p.v

     mdeg p == 
       p case R => 0$N
       degree(p.ts)$D

     init p == 
       p case R => error" Error in init from NSMP : #1 has no variables."
       leadingCoefficient(p.ts)$D

     head p == 
       p case R => p
       ([p.v,leadingMonomial(p.ts)$D]$VPoly)::Rep

     tail p == 
       p case R => 0$$
       red := reductum(p.ts)$D
       ground?(red)$D => (ground(red)$D)::Rep
       ([p.v,red]$VPoly)::Rep

     iteratedInitials p == 
       p case R => [] 
       p := leadingCoefficient(p.ts)$D
       cons(p,iteratedInitials(p)) 

     localDeepestInitial (p : $) : $ == 
       p case R => p 
       localDeepestInitial leadingCoefficient(p.ts)$D 

     deepestInitial p == 
       p case R => 
         error"Error in deepestInitial from NSMP : #1 has no variables."
       localDeepestInitial leadingCoefficient(p.ts)$D  

     mainMonomial p == 
       zero? p => 
         error"Error in mainMonomial from NSMP : the argument is zero"
       p case R => 1$$ 
       monomial(1$$,p.v,degree(p.ts)$D)

     leastMonomial p == 
       zero? p => 
         error"Error in leastMonomial from NSMP : the argument is zero"
       p case R => 1$$
       monomial(1$$,p.v,minimumDegree(p.ts)$D)

     mainCoefficients p == 
       zero? p => 
         error"Error in mainCoefficients from NSMP : the argument is zero"
       p case R => [p]
       coefficients(p.ts)$D

     leadingCoefficient(p:$,x:VarSet):$ == 
       (p case R) => p
       p.v = x => leadingCoefficient(p.ts)$D
       zero? (d := degree(p,x)) => p
       coefficient(p,x,d)

     localMonicModulo(a:$,b:$):$ ==
       -- b is assumed to have initial 1
       a case R => a
       a.v < b.v => a 
       mM: $
       if a.v > b.v
         then 
           m : D := map((a1:%):% +-> localMonicModulo(a1,b),a.ts)$SUP2
         else
           m : D := monicModulo(a.ts,b.ts)$D
       if ground?(m)$D 
          then 
            mM := (ground(m)$D)::Rep 
          else 
            mM := ([a.v,m]$VPoly)::Rep
       mM

     monicModulo (a,b) == 
       b case R => error"Error in monicModulo from NSMP : #2 is constant"
       ib : $ := init(b)@$
       not ground?(ib)$$ => 
         error"Error in monicModulo from NSMP : #2 is not monic"
       mM : $
       if not ((ib) = 1)$$
         then
           r : R := ground(ib)$$
           rec : Union(R,"failed"):= recip(r)$R
           (rec case "failed") =>
             error"Error in monicModulo from NSMP : #2 is not monic"
           a case R => a
           a := (rec::R) * a
           b := (rec::R) * b
           mM := ib * localMonicModulo (a,b)
         else
           mM := localMonicModulo (a,b)
       mM

     prem(a:$, b:$): $ == 
       -- with pseudoRemainder$NSUP
       b case R =>
         error "in prem$NSMP: ground? #2"
       db: N := degree(b.ts)$D
       lcb: $ := leadingCoefficient(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       delta: Z := max(test + 1$Z, 0$Z)
       (a case R) or (a.v < b.v) => lcb ** (delta::N) * a
       a.v = b.v =>
         r: D := pseudoRemainder(a.ts,b.ts)$D
         ground?(r) => return (ground(r)$D)::Rep 
         ([a.v,r]$VPoly)::Rep
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         delta := delta - 1$Z 
         test := degree(a,b.v)::Z - db
       lcb ** (delta::N) * a

     pquo (a:$, b:$)  : $ == 
       cPS := lazyPseudoDivide (a,b)
       c := (cPS.coef) ** (cPS.gap)
       c * cPS.quotient

     pseudoDivide(a:$, b:$): Record (quotient : $, remainder : $) ==
       -- from RPOLCAT
       cPS := lazyPseudoDivide(a,b)
       c := (cPS.coef) ** (cPS.gap)
       [c * cPS.quotient, c * cPS.remainder]

     lazyPrem(a:$, b:$): $ == 
       -- with lazyPseudoRemainder$NSUP
       -- Uses leadingCoefficient: ($, V) -> $
       b case R =>
         error "in lazyPrem$NSMP: ground? #2"
       (a case R) or (a.v < b.v) =>  a
       a.v = b.v => PSimp(lazyPseudoRemainder(a.ts,b.ts)$D,a.v)
       db: N := degree(b.ts)$D
       lcb: $ := leadingCoefficient(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         test := degree(a,b.v)::Z - db
       a

     lazyPquo (a:$, b:$) : $ ==
       -- with lazyPseudoQuotient$NSUP
       b case R =>
         error " in lazyPquo$NSMP: #2 is conctant"
       (a case R) or (a.v < b.v) => 0
       a.v = b.v => PSimp(lazyPseudoQuotient(a.ts,b.ts)$D,a.v)
       db: N := degree(b.ts)$D
       lcb: $ := leadingCoefficient(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       q := 0$$
       test: Z := degree(a,b.v)::Z - db
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         q := lcb * q + term
         test := degree(a,b.v)::Z - db
       q

     lazyPseudoDivide(a:$, b:$): _
         Record(coef:$, gap: N,quotient:$, remainder:$) == 
       -- with lazyPseudoDivide$NSUP
       b case R =>
         error " in lazyPseudoDivide$NSMP: #2 is conctant"
       (a case R) or (a.v < b.v) => [1$$,0$N,0$$,a]
       a.v = b.v =>
         cgqr := lazyPseudoDivide(a.ts,b.ts)
         [cgqr.coef, cgqr.gap, PSimp(cgqr.quotient,a.v), _
          PSimp(cgqr.remainder,a.v)]
       db: N := degree(b.ts)$D
       lcb: $ := leadingCoefficient(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       q := 0$$
       delta: Z := max(test + 1$Z, 0$Z) 
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         q := lcb * q + term
         delta := delta - 1$Z 
         test := degree(a,b.v)::Z - db
       [lcb, (delta::N), q, a]

     lazyResidueClass(a:$, b:$): Record(polnum:$, polden:$, power:N) == 
       -- with lazyResidueClass$NSUP
       b case R =>
         error " in lazyResidueClass$NSMP: #2 is conctant"
       lcb: $ := leadingCoefficient(b.ts)$D
       (a case R) or (a.v < b.v) => [a,lcb,0]
       a.v = b.v =>
         lrc := lazyResidueClass(a.ts,b.ts)$D
         [PSimp(lrc.polnum,a.v), lrc.polden, lrc.power]
       db: N := degree(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       pow: N := 0
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         pow := pow + 1
         test := degree(a,b.v)::Z - db
       [a, lcb, pow]

     if R has IntegralDomain
     then

       packD := PseudoRemainderSequence($,D)

       exactQuo(x:$, y:$):$ == 
         ex: Union($,"failed") := x exquo$$ y
         (ex case $) => ex::$
         error "in exactQuotient$NSMP: bad args"

       LazardQuotient(x:$, y:$, n: N):$ == 
         zero?(n) => error("LazardQuotient$NSMP : n = 0")
         (n = 1) => x
         a: N := 1
         while n >= (b := 2*a) repeat a := b
         c: $ := x
         n := (n - a)::N
         repeat       
           (a = 1) => return c
           a := a quo 2
           c := exactQuo(c*c,y)
           if n >= a then ( c := exactQuo(c*x,y) ; n := (n - a)::N )

       LazardQuotient2(p:$, a:$, b:$, n: N) ==
         zero?(n) => error " in LazardQuotient2$NSMP: bad #4"
         (n = 1) => p
         c: $  := LazardQuotient(a,b,(n-1)::N)
         exactQuo(c*p,b)

       next_subResultant2(p:$, q:$, z:$, s:$) ==
         PSimp(next_sousResultant2(p.ts,q.ts,z.ts,s)$packD,p.v)

       subResultantGcd(a:$, b:$): $ ==
         (a case R) or (b case R) => 
           error "subResultantGcd$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "subResultantGcd$NSMP: mvar(#1) ~= mvar(#2)"
         PSimp(subResultantGcd(a.ts,b.ts),a.v)

       halfExtendedSubResultantGcd1(a:$,b:$): Record (gcd: $, coef1: $) ==
         (a case R) or (b case R) => 
           error "halfExtendedSubResultantGcd1$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "halfExtendedSubResultantGcd1$NSMP: mvar(#1) ~= mvar(#2)"
         hesrg := halfExtendedSubResultantGcd1(a.ts,b.ts)$D
         [PSimp(hesrg.gcd,a.v), PSimp(hesrg.coef1,a.v)]

       halfExtendedSubResultantGcd2(a:$,b:$): Record (gcd: $, coef2: $) ==
         (a case R) or (b case R) => 
           error "halfExtendedSubResultantGcd2$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "halfExtendedSubResultantGcd2$NSMP: mvar(#1) ~= mvar(#2)"
         hesrg := halfExtendedSubResultantGcd2(a.ts,b.ts)$D
         [PSimp(hesrg.gcd,a.v), PSimp(hesrg.coef2,a.v)]

       extendedSubResultantGcd(a:$,b:$): Record (gcd: $, coef1: $, coef2: $) ==
         (a case R) or (b case R) => 
           error "extendedSubResultantGcd$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "extendedSubResultantGcd$NSMP: mvar(#1) ~= mvar(#2)"
         esrg := extendedSubResultantGcd(a.ts,b.ts)$D
         [PSimp(esrg.gcd,a.v),PSimp(esrg.coef1,a.v),PSimp(esrg.coef2,a.v)]  

       resultant(a:$, b:$): $ ==
         (a case R) or (b case R) => 
           error "resultant$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "resultant$NSMP: mvar(#1) ~= mvar(#2)"
         resultant(a.ts,b.ts)$D

       subResultantChain(a:$, b:$): List $ ==
         (a case R) or (b case R) => 
           error "subResultantChain$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "subResultantChain$NSMP: mvar(#1) ~= mvar(#2)"
         [PSimp(up,a.v) for up in subResultantsChain(a.ts,b.ts)]

       lastSubResultant(a:$, b:$): $ ==
         (a case R) or (b case R) => 
           error "lastSubResultant$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "lastSubResultant$NSMP: mvar(#1) ~= mvar(#2)"
         PSimp(lastSubResultant(a.ts,b.ts),a.v)

       if R has EuclideanDomain
       then

         exactQuotient (a:$,b:R) ==
           (b = 1) => a
           a case R => (a::R quo$R b)::$
           ([a.v, map((a1:%):% +-> exactQuotient(a1,b),a.ts)$SUP2]$VPoly)::Rep

         exactQuotient! (a:$,b:R) ==
           (b = 1) => a
           a case R => (a::R quo$R b)::$
           a.ts := map((a1:%):% +-> exactQuotient!(a1,b),a.ts)$SUP2
           a

       else

         exactQuotient (a:$,b:R) ==
           (b = 1) => a
           a case R => ((a::R exquo$R b)::R)::$
           ([a.v, map((a1:%):% +-> exactQuotient(a1,b),a.ts)$SUP2]$VPoly)::Rep

         exactQuotient! (a:$,b:R) == 
           (b = 1) => a
           a case R => ((a::R exquo$R b)::R)::$
           a.ts := map((a1:%):% +-> exactQuotient!(a1,b),a.ts)$SUP2
           a

     if R has GcdDomain
     then

       localGcd(r:R,p:$):R ==
         p case R => gcd(r,p::R)$R
         gcd(r,content(p))$R         

       gcd(r:R,p:$):R ==
         (r = 1) => r
         zero? p => r
         localGcd(r,p)

       content p ==
         p case R => p
         up : D := p.ts
         r := 0$R
         while (not zero? up) and (not (r = 1)) repeat
           r := localGcd(r,leadingCoefficient(up))
           up := reductum up
         r

       primitivePart! p ==
         zero? p => p
         p case R => 1$$
         cp := content(p)
         p.ts := 
           unitCanonical(map((a1:%):% +-> exactQuotient!(a1,cp),p.ts)$SUP2)$D
         p