Edit detail for SandBox Polynomials revision 1 of 2

changed:
-
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.

\begin{spad}
)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))
      differentiate(p: %,v: OV) ==
         R has PDRING SYMBOL => multivariate(differentiate(univariate(p),convert(v)@Symbol)$SUP(R),v)
         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(#1.k > #2.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(#1 > #2,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)

\end{spad}

Test:
\begin{axiom}
P := DMP([x,y], EXPR INT) 
a :P := x 
b := a/x
univariate(b)
differentiate(b,x)
\end{axiom}

From wyscc Sat Feb 25 23:14:40 -0600 2006
From: wyscc
Date: Sat, 25 Feb 2006 23:14:40 -0600
Subject: Quizzes
Message-ID: <20060225231440-0600@wiki.axiom-developer.org>


Instructions: What is the output for each of the following?
Do the problem only one at a time. You should write down your answer, its type, or possible errors, and with explanations. Then check it (using )set mess bot off), revise your explanation if necessary, or explain Axiom's output. Finally, use )set mess bot on to see exactly what the Interpreter did. Then and only then, advance to the next problem.
<pre>
  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)
</pre>
Answer: see SandBoxPolynomialQuiz

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.

spad

)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))
      differentiate(p: %,v: OV) ==
         R has PDRING SYMBOL => multivariate(differentiate(univariate(p),convert(v)@Symbol)$SUP(R),v)
         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(#1.k > #2.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(#1 > #2,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)

spad

   Compiling FriCAS source code from file 
      /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/6733095636270387759-25px001.spad
      using old system compiler.
   GDMP abbreviates domain GeneralDistributedMultivariatePolynomial 
******** Spad syntax error detected ********
The prior line was:

  139>               maxdeg.i := max(maxdeg.i, tdeg.i)

The current line is:

  140>          [index(i:PositiveInteger) for i in 1..n | maxdeg.i^=0]

The number of valid tokens is 2.
The prior token was #S(TOKEN :SYMBOL ^ :TYPE KEYWORD :NONBLANK 140)
The current token is #S(TOKEN :SYMBOL = :TYPE KEYWORD :NONBLANK 140)
The next token is #S(TOKEN :SYMBOL 0 :TYPE NUMBER :NONBLANK 140)

Test:

fricas

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

(1)

Type: Type

fricas

a :P := x

(2)

Type: DistributedMultivariatePolynomial?([x,y],Expression(Integer))

fricas

b := a/x

(3)

Type: DistributedMultivariatePolynomial?([x,y],Expression(Integer))

fricas

univariate(b)

(4)

Type: SparseUnivariatePolynomial?(Expression(Integer))

fricas

differentiate(b,x)

(5)

Type: DistributedMultivariatePolynomial?([x,y],Expression(Integer))

Quizzes --wyscc, Sat, 25 Feb 2006 23:14:40 -0600 reply

Instructions: What is the output for each of the following? Do the problem only one at a time. You should write down your answer, its type, or possible errors, and with explanations. Then check it (using )set mess bot off), revise your explanation if necessary, or explain Axiom's output. Finally, use )set mess bot on to see exactly what the Interpreter did. Then and only then, advance to the next problem.

  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)

Answer: see SandBoxPolynomialQuiz?