Testing a change from:

      differentiate(p: %,v: OV) ==
            multivariate(differentiate(univariate(p,v)),v)

to:

      differentiate(p: %,v: OV) ==
            multivariate(differentiate(univariate(p,v),v),v)

i.e. differentiate should be extended from the underlying Ring.

)abbrev domain GDMP GeneralDistributedMultivariatePolynomial
++ Author: Barry Trager
++ Date Created:
++ Date Last Updated:
++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
++ resultant, gcd, leadingCoefficient
++ Related Constructors: DistributedMultivariatePolynomial,
++ HomogeneousDistributedMultivariatePolynomial
++ Also See: Polynomial
++ AMS Classifications:
++ Keywords: polynomial, multivariate, distributed
++ References:
++ Description:
++   This type supports distributed multivariate polynomials
++ whose variables are from a user specified list of symbols.
++ The coefficient ring may be non commutative,
++ but the variables are assumed to commute.
++ The term ordering is specified by its third parameter.
++ Suggested types which define term orderings include: \spadtype{DirectProduct},
++ \spadtype{HomogeneousDirectProduct}, \spadtype{SplitHomogeneousDirectProduct}
++ and finally \spadtype{OrderedDirectProduct} which accepts an arbitrary user
++ function to define a term ordering.

GeneralDistributedMultivariatePolynomial(vl,R,E): public == private where
  vl: List Symbol
  R: Ring
  E: DirectProductCategory(#vl,NonNegativeInteger)
  OV  ==> OrderedVariableList(vl)
  SUP ==> SparseUnivariatePolynomial
  NNI ==> NonNegativeInteger

  public == PolynomialCategory(R,E,OV) with
      reorder: (%,List Integer) -> %
        ++ reorder(p, perm) applies the permutation perm to the variables
        ++ in a polynomial and returns the new correctly ordered polynomial
      univariate:(%,OV) -> SUP(R)

  private == PolynomialRing(R,E) add
    --representations
      Term := Record(k:E,c:R)
      Rep := List Term
      n := #vl
      Vec ==> Vector(NonNegativeInteger)
      zero?(p : %): Boolean == null(p::Rep)

      totalDegree p ==
         zero? p => 0
         "max"/[reduce("+",(t.k)::(Vector NNI), 0) for t in p]

      monomial(p:%, v: OV,e: NonNegativeInteger):% ==
         locv := lookup v
         p*monomial(1,
            directProduct [if z=locv then e else 0 for z in 1..n]$Vec)

      coerce(v: OV):% == monomial(1,v,1)

      listCoef(p : %): List R ==
        rec : Term
        [rec.c for rec in (p::Rep)]

      mainVariable(p: %) ==
         zero?(p) => "failed"
         for v in vl repeat
           vv := variable(v)::OV
           if degree(p,vv)>0 then return vv
         "failed"

      ground?(p) == mainVariable(p) case "failed"

      retract(p : %): R ==
          not ground? p => error "not a constant"
          leadingCoefficient p

      retractIfCan(p : %): Union(R,"failed") ==
        ground?(p) => leadingCoefficient p
        "failed"

      degree(p: %,v: OV) == degree(univariate(p,v))
      minimumDegree(p: %,v: OV) == minimumDegree(univariate(p,v))
      if R has PartialDifferentialRing(Symbol) then

         differentiate(p : %, v : OV) ==
            pu := univariate(p)
            dpu := map((x : R) : R +-> differentiate(x, convert(v)@Symbol), pu)
            multivariate(differentiate(pu) + dpu, v)
      else

         differentiate(p : %, v : OV) ==
             multivariate(differentiate(univariate(p,v)),v)

      degree(p: %,lv: List OV) == [degree(p,v) for v in lv]
      minimumDegree(p: %,lv: List OV) == [minimumDegree(p,v) for v in lv]

      numberOfMonomials(p:%) ==
        l : Rep := p::Rep
        null(l) => 1
        #l

      monomial?(p : %): Boolean ==
        l : Rep := p::Rep
        null(l) or null rest(l)

      if R has OrderedRing then
        maxNorm(p : %): R ==
          l : List R := nil
          r,m : R
          m := 0
          for r in listCoef(p) repeat
            if r > m then m := r
            else if (-r) > m then m := -r
          m

      --trailingCoef(p : %) ==
      --  l : Rep := p : Rep
      --  null l => 0

      --  r : Term := last l
      --  r.c

      --leadingPrimitiveMonomial(p : %) ==
      --  ground?(p) => 1$%
      --  r : Term := first(p:Rep)
      --  r := [r.k,1$R]$Term     -- new cell

      -- list(r)$Rep :: %

    -- The following 2 defs are inherited from PolynomialRing

      --leadingMonomial(p : %) ==
      --  ground?(p) => p
      --  r : Term := first(p:Rep)
      --  r := [r.k,r.c]$Term     -- new cell
      --  list(r)$Rep :: %

      --reductum(p : %): % ==
      --  ground? p => 0$%
      --  (rest(p:Rep)):%

      if R has Field then
        (p : %) / (r : R) == inv(r) * p

      variables(p: %) ==
         maxdeg:Vector(NonNegativeInteger) := new(n,0)
         while not zero?(p) repeat
            tdeg := degree p
            p := reductum p
            for i in 1..n repeat
              maxdeg.i := max(maxdeg.i, tdeg.i)
         [index(i::PositiveInteger) for i in 1..n | maxdeg.i ~= 0]

      reorder(p: %,perm: List Integer):% ==
         #perm ~= n => error "must be a complete permutation of all vars"
         q := [[directProduct [term.k.j for j in perm]$Vec,term.c]$Term
                         for term in p]
         sort((x, y) +-> x.k > x.k, q)

      --coerce(dp:DistributedMultivariatePolynomial(vl,R)):% ==
      --   q:=dp:List(Term)
      --   sort(#1.k > #2.k,q):%

      univariate(p: %,v: OV):SUP(%) ==
         zero?(p) => 0
         exp := degree p
         locv := lookup v
         deg:NonNegativeInteger := 0
         nexp := directProduct [if i=locv then (deg :=exp.i;0) else exp.i
                                        for i in 1..n]$Vec
         monomial(monomial(leadingCoefficient p,nexp),deg)+
                      univariate(reductum p,v)

      eval(p: %,v: OV,val:%):% == univariate(p,v)(val)

      eval(p: %,v: OV,val:R):% == eval(p,v,val::%)$%

      eval(p: %,lv: List OV,lval: List R):% ==
         lv = [] => p
         eval(eval(p,first lv,(first lval)::%)$%, rest lv, rest lval)$%

      -- assume Lvar are sorted correctly
      evalSortedVarlist(p: %,Lvar: List OV,Lpval: List %):% ==
        v := mainVariable p
        v case "failed" => p
        pv := v:: OV
        Lvar=[] or Lpval=[] => p
        mvar := Lvar.first
        mvar > pv => evalSortedVarlist(p,Lvar.rest,Lpval.rest)
        pval := Lpval.first
        pts:SUP(%):= map(evalSortedVarlist(#1,Lvar,Lpval),univariate(p,pv))
        mvar=pv => pts(pval)
        multivariate(pts,pv)

      eval(p:%,Lvar:List OV,Lpval:List %) ==
        nlvar:List OV := sort((x : OV, y : OV) : Boolean +-> x > y, Lvar)
        nlpval :=
           Lvar = nlvar => Lpval
           nlpval := [Lpval.position(mvar,Lvar) for mvar in nlvar]
        evalSortedVarlist(p,nlvar,nlpval)

      multivariate(p1:SUP(%),v: OV):% ==
        0=p1 => 0
        degree p1 = 0 => leadingCoefficient p1
        leadingCoefficient(p1)*(v::%)^degree(p1) +
                  multivariate(reductum p1,v)

      univariate(p: %):SUP(R) ==
        (v := mainVariable p) case "failed" =>
                      monomial(leadingCoefficient p,0)
        q := univariate(p,v:: OV)
        ans:SUP(R) := 0
        while q ~= 0 repeat
          ans := ans + monomial(ground leadingCoefficient q,degree q)
          q := reductum q
        ans

      multivariate(p:SUP(R),v: OV):% ==
        0=p => 0
        (leadingCoefficient p)*monomial(1,v,degree p) +
                       multivariate(reductum p,v)

      if R has GcdDomain then
        content(p: %):R ==
          zero?(p) => 0
          "gcd"/[t.c for t in p]

        if R has EuclideanDomain and not(R has FloatingPointSystem)  then
          gcd(p: %,q:%):% ==
            gcd(p,q)$PolynomialGcdPackage(E,OV,R,%)

        else gcd(p: %,q:%):% ==
            r : R
            (pv := mainVariable(p)) case "failed" =>
              (r := leadingCoefficient p) = 0$R => q
              gcd(r,content q)::%
            (qv := mainVariable(q)) case "failed" =>
              (r := leadingCoefficient q) = 0$R => p
              gcd(r,content p)::%
            pv<qv => gcd(p,content univariate(q,qv))
            qv<pv => gcd(q,content univariate(p,pv))
            multivariate(gcd(univariate(p,pv),univariate(q,qv)),pv)

      coerce(p: %) : OutputForm ==
        zero?(p) => (0$R) :: OutputForm
        l,lt : List OutputForm
        lt := nil
        vl1 := [v::OutputForm for v in vl]
        for t in reverse p repeat
          l := nil
          for i in 1..#vl1 repeat
            t.k.i = 0 => l
            t.k.i = 1 => l := cons(vl1.i,l)
            l := cons(vl1.i ^ t.k.i ::OutputForm,l)
          l := reverse l
          if (t.c ~= 1) or (null l) then l := cons(t.c :: OutputForm,l)
          1 = #l => lt := cons(first l,lt)
          lt := cons(reduce("*",l),lt)
        1 = #lt => first lt
        reduce("+",lt)

   Compiling FriCAS source code from file 
      /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/6560064616104104256-25px001.spad
      using old system compiler.
   GDMP abbreviates domain GeneralDistributedMultivariatePolynomial 
------------------------------------------------------------------------
   initializing NRLIB GDMP for GeneralDistributedMultivariatePolynomial 
   compiling into NRLIB GDMP 
   processing macro definition Vec ==> Vector NonNegativeInteger 
   compiling exported zero? : $ -> Boolean
****** comp fails at level 1 with expression: ******
error in function zero? 

((|null| (|::| |p| |Rep|)))
****** level 1  ******
$x:= (null (:: p Rep))
$m:= (Boolean)
$f:=
((((|p| # #) (|Vec| #) (|n| #) (* #) ...)))

   >> Apparent user error:
   cannot compile (null (:: p Rep))

Test:

P := DMP([x,y], EXPR INT)

(1)

a :P := x

(2)

b := a/x

(3)

univariate(b)

(4)

differentiate(b,x)

(5)

  variables (2*x+1/x)$DMP([x], EXPR INT)
  variables (2*y+1/y)$DMP([y], INT)
  a:=(2*x+1/x)$DMP([x], EXPR INT); variables a
  b:=(2*y+1/y)$DMP([y], INT); variables b
  x:DMP([x], EXPR INT):=x; variables (2*x+1/x)
  y:DMP([y], INT):=y; variables (2*y+1/y)