axiom-developer
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Axiom-developer] 20081211.01.tpd.patch (regression test suite cleanup)


From: daly
Subject: [Axiom-developer] 20081211.01.tpd.patch (regression test suite cleanup)
Date: Fri, 12 Dec 2008 13:32:31 -0600

The patch cleans up some minor breakage in the regression test suite.

It also fixes cwmmt in preparation for removing noweb from the tool chain.

Tim

=====================================================================
diff --git a/books/bookvol10.3.pamphlet b/books/bookvol10.3.pamphlet
index 3510eef..33e4bb1 100644
--- a/books/bookvol10.3.pamphlet
+++ b/books/bookvol10.3.pamphlet
@@ -15348,6 +15348,117 @@ Equation(S: Type): public == private where
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain EMR EuclideanModularRing}
+\pagehead{EuclideanModularRing}{EMR}
+\pagepic{ps/v103euclideanmodularring.ps}{EMR}{1.00}
+See also:\\
+\refto{ModularRing}{MODRING}
+\refto{ModularField}{MODFIELD}
+<<domain EMR EuclideanModularRing>>=
+)abbrev domain EMR EuclideanModularRing
+++ Description:
+++ These domains are used for the factorization and gcds
+++ of univariate polynomials over the integers in order to work modulo
+++ different  primes.
+++ See \spadtype{ModularRing}, \spadtype{ModularField}
+EuclideanModularRing(S,R,Mod,reduction:(R,Mod) -> R,
+                     merge:(Mod,Mod) -> Union(Mod,"failed"),
+                      exactQuo : (R,R,Mod) -> Union(R,"failed")) : C == T
+ where
+  S    :  CommutativeRing
+  R    :  UnivariatePolynomialCategory S
+  Mod  :  AbelianMonoid
+
+  C == EuclideanDomain with
+                modulus :   %     -> Mod
+                       ++ modulus(x) \undocumented
+                coerce  :   %     -> R
+                       ++ coerce(x) \undocumented
+                reduce  : (R,Mod) -> %
+                       ++ reduce(r,m) \undocumented
+                exQuo   :  (%,%)  -> Union(%,"failed")
+                       ++ exQuo(x,y) \undocumented
+                recip   :    %    -> Union(%,"failed")
+                       ++ recip(x) \undocumented
+                inv     :    %    -> %
+                       ++ inv(x) \undocumented
+                elt     : (%, R)  -> R
+                       ++ elt(x,r) or x.r \undocumented
+
+  T == ModularRing(R,Mod,reduction,merge,exactQuo) add
+
+    --representation
+      Rep:= Record(val:R,modulo:Mod)
+    --declarations
+      x,y,z: %
+
+      divide(x,y) ==
+        t:=merge(x.modulo,y.modulo)
+        t case "failed" => error "incompatible moduli"
+        xm:=t::Mod
+        yv:=y.val
+        invlcy:R
+--        if one? leadingCoefficient yv then invlcy:=1
+        if (leadingCoefficient yv = 1) then invlcy:=1
+        else
+          invlcy:=(inv reduce((leadingCoefficient yv)::R,xm)).val
+          yv:=reduction(invlcy*yv,xm)
+        r:=monicDivide(x.val,yv)
+        [reduce(invlcy*r.quotient,xm),reduce(r.remainder,xm)]
+
+      if R has fmecg:(R,NonNegativeInteger,S,R)->R
+         then x rem y  ==
+           t:=merge(x.modulo,y.modulo)
+           t case "failed" => error "incompatible moduli"
+           xm:=t::Mod
+           yv:=y.val
+           invlcy:R
+--           if not one? leadingCoefficient yv then
+           if not (leadingCoefficient yv = 1) then
+             invlcy:=(inv reduce((leadingCoefficient yv)::R,xm)).val
+             yv:=reduction(invlcy*yv,xm)
+           dy:=degree yv
+           xv:=x.val
+           while (d:=degree xv - dy)>=0 repeat
+                 xv:=reduction(fmecg(xv,d::NonNegativeInteger,
+                                     leadingCoefficient xv,yv),xm)
+                 xv = 0 => return [xv,xm]$Rep
+           [xv,xm]$Rep
+         else x rem y  == 
+           t:=merge(x.modulo,y.modulo)
+           t case "failed" => error "incompatible moduli"
+           xm:=t::Mod
+           yv:=y.val
+           invlcy:R
+--           if not one? leadingCoefficient yv then
+           if not (leadingCoefficient yv = 1) then
+             invlcy:=(inv reduce((leadingCoefficient yv)::R,xm)).val
+             yv:=reduction(invlcy*yv,xm)
+           r:=monicDivide(x.val,yv)
+           reduce(r.remainder,xm)
+
+      euclideanSize x == degree x.val
+
+      unitCanonical x ==
+        zero? x => x
+        degree(x.val) = 0 => 1
+--        one? leadingCoefficient(x.val) => x
+        (leadingCoefficient(x.val) = 1) => x
+        invlcx:%:=inv reduce((leadingCoefficient(x.val))::R,x.modulo)
+        invlcx * x
+
+      unitNormal x ==
+--        zero?(x) or one?(leadingCoefficient(x.val)) => [1, x, 1]
+        zero?(x) or ((leadingCoefficient(x.val)) = 1) => [1, x, 1]
+        lcx := reduce((leadingCoefficient(x.val))::R,x.modulo)
+        invlcx:=inv lcx
+        degree(x.val) = 0 => [lcx, 1, invlcx]
+        [lcx, invlcx * x, invlcx]
+
+      elt(x : %,s : R) : R == reduction(elt(x.val,s),x.modulo)
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain EXPEXPAN ExponentialExpansion}
 \pagehead{ExponentialExpansion}{EXPEXPAN}
 \pagepic{ps/v103exponentialexpansion.ps}{EXPEXPAN}{1.00}
@@ -26690,6 +26801,82 @@ GeneralDistributedMultivariatePolynomial(vl,R,E): 
public == private where
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain GMODPOL GeneralModulePolynomial}
+\pagehead{GeneralModulePolynomial}{GMODPOL}
+\pagepic{ps/v103generalmodulepolynomial.ps}{GMODPOL}{1.00}
+See also:\\
+\refto{ModuleMonomial}{MODMONOM}
+<<domain GMODPOL GeneralModulePolynomial>>=
+)abbrev domain GMODPOL GeneralModulePolynomial
+++ Description:
+++ This package \undocumented
+GeneralModulePolynomial(vl, R, IS, E, ff, P): public  ==  private where
+  vl: List(Symbol)
+  R: CommutativeRing
+  IS: OrderedSet
+  NNI ==> NonNegativeInteger
+  E: DirectProductCategory(#vl, NNI)
+  MM ==> Record(index:IS, exponent:E)
+  ff: (MM, MM) -> Boolean
+  OV  ==> OrderedVariableList(vl)
+  P: PolynomialCategory(R, E, OV)
+  ModMonom ==> ModuleMonomial(IS, E, ff)
+
+
+  public  ==  Join(Module(P), Module(R))  with
+        leadingCoefficient: $ -> R
+               ++ leadingCoefficient(x) \undocumented
+        leadingMonomial: $ -> ModMonom
+               ++ leadingMonomial(x) \undocumented
+        leadingExponent: $ -> E
+               ++ leadingExponent(x) \undocumented
+        leadingIndex: $ -> IS
+               ++ leadingIndex(x) \undocumented
+        reductum: $ -> $
+               ++ reductum(x) \undocumented
+        monomial: (R, ModMonom) -> $
+               ++ monomial(r,x) \undocumented
+        unitVector: IS -> $
+               ++ unitVector(x) \undocumented
+        build: (R, IS, E) -> $
+               ++ build(r,i,e) \undocumented
+        multMonom: (R, E, $) -> $
+               ++ multMonom(r,e,x) \undocumented
+        "*": (P,$) -> $
+               ++ p*x \undocumented
+
+
+  private  ==  FreeModule(R, ModMonom)  add
+        Rep:= FreeModule(R, ModMonom)
+        leadingMonomial(p:$):ModMonom == leadingSupport(p)$Rep
+        leadingExponent(p:$):E == exponent(leadingMonomial p)
+        leadingIndex(p:$):IS == index(leadingMonomial p)
+        unitVector(i:IS):$ == monomial(1,[i, 0$E]$ModMonom)
+
+
+ -----------------------------------------------------------------------------
+
+        build(c:R, i:IS, e:E):$  ==  monomial(c, construct(i, e))
+
+ -----------------------------------------------------------------------------
+
+     ----   WARNING: assumes c ^= 0
+
+        multMonom(c:R, e:E, mp:$):$  ==
+            zero? mp => mp
+            monomial(c * leadingCoefficient mp, [leadingIndex mp,
+                     e + leadingExponent mp]) + multMonom(c, e, reductum mp)
+
+ -----------------------------------------------------------------------------
+
+
+        ((p:P) * (mp:$)):$  ==
+            zero? p => 0
+            multMonom(leadingCoefficient p, degree p, mp) +
+               reductum(p) * mp
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain GCNAALG GenericNonAssociativeAlgebra}
 \pagehead{GenericNonAssociativeAlgebra}{GCNAALG}
 \pagepic{ps/v103genericnonassociativealgebra.ps}{GCNAALG}{1.00}
@@ -28243,6 +28430,47 @@ 
IndexedDirectProductOrderedAbelianMonoidSup(A:OrderedAbelianMonoidSup,S:OrderedS
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain INDE IndexedExponents}
+\pagehead{IndexedExponents}{INDE}
+\pagepic{ps/v103indexedexponents.ps}{INDE}{1.00}
+See also:\\
+\refto{Polynomial}{POLY}
+\refto{MultivariatePolynomial}{MPOLY}
+\refto{SparseMultivariatePolynomial}{SMP}
+<<domain INDE IndexedExponents>>=
+)abbrev domain INDE IndexedExponents
+++ Author: James Davenport
+++ Date Created:
+++ Date Last Updated:
+++ Basic Functions:
+++ Related Constructors:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ References:
+++ Description:
+++   IndexedExponents of an ordered set of variables gives a representation
+++ for the degree of polynomials in commuting variables. It gives an ordered
+++ pairing of non negative integer exponents with variables
+
+IndexedExponents(Varset:OrderedSet): C == T where
+  C == Join(OrderedAbelianMonoidSup,
+            IndexedDirectProductCategory(NonNegativeInteger,Varset))
+  T == IndexedDirectProductOrderedAbelianMonoidSup(NonNegativeInteger,Varset) 
add
+      Term:=  Record(k:Varset,c:NonNegativeInteger)
+      Rep:=  List Term
+      x:%
+      t:Term
+      coerceOF(t):OutputForm ==     --++ converts term to OutputForm
+         t.c = 1 => (t.k)::OutputForm
+         (t.k)::OutputForm ** (t.c)::OutputForm
+      coerce(x):OutputForm == ++ converts entire exponents to OutputForm
+         null x => 1::Integer::OutputForm
+         null rest x => coerceOF(first x)
+         reduce("*",[coerceOF t for t in x])
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain IFARRAY IndexedFlexibleArray}
 <<dot>>=
 "IFARRAY" -> "A1AGG"
@@ -28694,6 +28922,81 @@ IndexedList(S:Type, mn:Integer): Exports == 
Implementation where
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain IMATRIX IndexedMatrix}
+\pagehead{IndexedMatrix}{IMATRIX}
+\pagepic{ps/v103indexedmatrix.ps}{IMATRIX}{1.00}
+See also:\\
+\refto{Matrix}{MATRIX}
+\refto{RectangularMatrix}{RMATRIX}
+\refto{SquareMatrix}{SQMATRIX}
+<<domain IMATRIX IndexedMatrix>>=
+)abbrev domain IMATRIX IndexedMatrix
+++ Author: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: Matrix, RectangularMatrix, SquareMatrix,
+++   StorageEfficientMatrixOperations
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   An \spad{IndexedMatrix} is a matrix where the minimal row and column
+++   indices are parameters of the type.  The domains Row and Col
+++   are both IndexedVectors.
+++   The index of the 'first' row may be obtained by calling the
+++   function \spadfun{minRowIndex}.  The index of the 'first' column may
+++   be obtained by calling the function \spadfun{minColIndex}.  The index of
+++   the first element of a 'Row' is the same as the index of the
+++   first column in a matrix and vice versa.
+IndexedMatrix(R,mnRow,mnCol): Exports == Implementation where
+  R : Ring
+  mnRow, mnCol : Integer
+  Row ==> IndexedVector(R,mnCol)
+  Col ==> IndexedVector(R,mnRow)
+  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,$)
+ 
+  Exports ==> MatrixCategory(R,Row,Col)
+ 
+  Implementation ==>
+    InnerIndexedTwoDimensionalArray(R,mnRow,mnCol,Row,Col) add
+ 
+      swapRows_!(x,i1,i2) ==
+        (i1 < minRowIndex(x)) or (i1 > maxRowIndex(x)) or _
+           (i2 < minRowIndex(x)) or (i2 > maxRowIndex(x)) =>
+             error "swapRows!: index out of range"
+        i1 = i2 => x
+        minRow := minRowIndex x
+        xx := x pretend PrimitiveArray(PrimitiveArray(R))
+        n1 := i1 - minRow; n2 := i2 - minRow
+        row1 := qelt(xx,n1)
+        qsetelt_!(xx,n1,qelt(xx,n2))
+        qsetelt_!(xx,n2,row1)
+        xx pretend $
+ 
+      if R has commutative("*") then
+ 
+        determinant x == determinant(x)$MATLIN
+        minordet    x == minordet(x)$MATLIN
+ 
+      if R has EuclideanDomain then
+ 
+        rowEchelon  x == rowEchelon(x)$MATLIN
+ 
+      if R has IntegralDomain then
+ 
+        rank        x == rank(x)$MATLIN
+        nullity     x == nullity(x)$MATLIN
+        nullSpace   x == nullSpace(x)$MATLIN
+ 
+      if R has Field then
+ 
+        inverse     x == inverse(x)$MATLIN
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain IARRAY1 IndexedOneDimensionalArray}
 <<dot>>=
 "IARRAY1" -> "A1AGG"
@@ -29225,6 +29528,198 @@ 
InnerIndexedTwoDimensionalArray(R,mnRow,mnCol,Row,Col):_
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain INFORM InputForm}
+\pagehead{InputForm}{INFORM}
+\pagepic{ps/v103inputform.ps}{INFORM}{1.00}
+<<domain INFORM InputForm>>=
+)abbrev domain INFORM InputForm
+++ Parser forms
+++ Author: Manuel Bronstein
+++ Date Created: ???
+++ Date Last Updated: 19 April 1991
+++ Description:
+++   Domain of parsed forms which can be passed to the interpreter.
+++   This is also the interface between algebra code and facilities
+++   in the interpreter.
+
+--)boot $noSubsumption := true
+
+InputForm():
+  Join(SExpressionCategory(String,Symbol,Integer,DoubleFloat,OutputForm),
+       ConvertibleTo SExpression) with
+    interpret: % -> Any
+      ++ interpret(f) passes f to the interpreter.
+    convert  : SExpression -> %
+      ++ convert(s) makes s into an input form.
+    binary   : (%, List %) -> %
+      ++ \spad{binary(op, [a1,...,an])} returns the input form
+      ++ corresponding to  \spad{a1 op a2 op ... op an}.
+      ++
+      ++X a:=[1,2,3]::List(InputForm)
+      ++X binary(_+::InputForm,a)
+
+    function : (%, List Symbol, Symbol) -> %
+      ++ \spad{function(code, [x1,...,xn], f)} returns the input form
+      ++ corresponding to \spad{f(x1,...,xn) == code}.
+    lambda   : (%, List Symbol) -> %
+      ++ \spad{lambda(code, [x1,...,xn])} returns the input form
+      ++ corresponding to \spad{(x1,...,xn) +-> code} if \spad{n > 1},
+      ++ or to \spad{x1 +-> code} if \spad{n = 1}.
+    "+"      : (%, %) -> %
+      ++ \spad{a + b} returns the input form corresponding to \spad{a + b}.
+    "*"      : (%, %) -> %
+      ++ \spad{a * b} returns the input form corresponding to \spad{a * b}.
+    "/"      : (%, %) -> %
+      ++ \spad{a / b} returns the input form corresponding to \spad{a / b}.
+    "**"     : (%, NonNegativeInteger) -> %
+      ++ \spad{a ** b} returns the input form corresponding to \spad{a ** b}.
+    "**"     : (%, Integer) -> %
+      ++ \spad{a ** b} returns the input form corresponding to \spad{a ** b}.
+    0        : constant -> %
+      ++ \spad{0} returns the input form corresponding to 0.
+    1        : constant -> %
+      ++ \spad{1} returns the input form corresponding to 1.
+    flatten  : % -> %
+      ++ flatten(s) returns an input form corresponding to s with
+      ++ all the nested operations flattened to triples using new
+      ++ local variables.
+      ++ If s is a piece of code, this speeds up
+      ++ the compilation tremendously later on.
+    unparse  : % -> String
+      ++ unparse(f) returns a string s such that the parser
+      ++ would transform s to f.
+      ++ Error: if f is not the parsed form of a string.
+    parse : String -> %
+      ++ parse is the inverse of unparse. It parses a string to InputForm.
+    declare  : List %   -> Symbol
+      ++ declare(t) returns a name f such that f has been
+      ++ declared to the interpreter to be of type t, but has
+      ++ not been assigned a value yet.
+      ++ Note: t should be created as \spad{devaluate(T)$Lisp} where T is the
+      ++ actual type of f (this hack is required for the case where
+      ++ T is a mapping type).
+    compile  : (Symbol, List %) -> Symbol
+      ++ \spad{compile(f, [t1,...,tn])} forces the interpreter to compile
+      ++ the function f with signature \spad{(t1,...,tn) -> ?}.
+      ++ returns the symbol f if successful.
+      ++ Error: if f was not defined beforehand in the interpreter,
+      ++ or if the ti's are not valid types, or if the compiler fails.
+ == SExpression add
+    Rep := SExpression
+
+    mkProperOp: Symbol -> %
+    strsym    : % -> String
+    tuplify   : List Symbol -> %
+    flatten0  : (%, Symbol, NonNegativeInteger) ->
+                                             Record(lst: List %, symb:%)
+
+    0                        == convert(0::Integer)
+    1                        == convert(1::Integer)
+    convert(x:%):SExpression == x pretend SExpression
+    convert(x:SExpression):% == x
+
+    conv(ll : List %): % ==
+      convert(ll pretend List SExpression)$SExpression pretend %
+
+    lambda(f,l) == conv([convert("+->"::Symbol),tuplify l,f]$List(%))
+
+    interpret x ==
+      v := interpret(x)$Lisp
+      mkObj(unwrap(objVal(v)$Lisp)$Lisp, objMode(v)$Lisp)$Lisp
+
+    convert(x:DoubleFloat):% ==
+      zero? x => 0
+--      one? x => 1
+      (x = 1) => 1
+      convert(x)$Rep
+
+    flatten s ==
+      -- will not compile if I use 'or'
+      atom? s => s
+      every?(atom?,destruct s)$List(%) => s
+      sy := new()$Symbol
+      n:NonNegativeInteger := 0
+      l2 := [flatten0(x, sy, n := n + 1) for x in rest(l := destruct s)]
+      conv(concat(convert("SEQ"::Symbol)@%,
+        concat(concat [u.lst for u in l2], conv(
+           [convert("exit"::Symbol)@%, 1$%, conv(concat(first l,
+               [u.symb for u in l2]))@%]$List(%))@%)))@%
+
+    flatten0(s, sy, n) ==
+      atom? s => [nil(), s]
+      a := convert(concat(string sy, convert(n)@String)::Symbol)@%
+      l2 := [flatten0(x, sy, n := n+1) for x in rest(l := destruct s)]
+      [concat(concat [u.lst for u in l2], conv([convert(
+        "LET"::Symbol)@%, a, conv(concat(first l,
+             [u.symb for u in l2]))@%]$List(%))@%), a]
+
+    strsym s ==
+      string? s => string s
+      symbol? s => string symbol s
+      error "strsym: form is neither a string or symbol"
+
+    unparse x ==
+      atom?(s:% := form2String(x)$Lisp) => strsym s
+      concat [strsym a for a in destruct s]
+
+    parse(s:String):% ==
+      ncParseFromString(s)$Lisp
+
+    declare signature ==
+      declare(name := new()$Symbol, signature)$Lisp
+      name
+
+    compile(name, types) ==
+      symbol car cdr car
+        selectLocalMms(mkProperOp name, convert(name)@%,
+          types, nil$List(%))$Lisp
+
+    mkProperOp name ==
+      op := mkAtree(nme := convert(name)@%)$Lisp
+      transferPropsToNode(nme, op)$Lisp
+      convert op
+
+    binary(op, args) ==
+      (n := #args) < 2 => error "Need at least 2 arguments"
+      n = 2 => convert([op, first args, last args]$List(%))
+      convert([op, first args, binary(op, rest args)]$List(%))
+
+    tuplify l ==
+      empty? rest l => convert first l
+      conv
+        concat(convert("Tuple"::Symbol), [convert x for x in l]$List(%))
+
+    function(f, l, name) ==
+      nn := convert(new(1 + #l, convert(nil()$List(%)))$List(%))@%
+      conv([convert("DEF"::Symbol), conv(cons(convert(name)@%,
+                        [convert(x)@% for x in l])), nn, nn, f]$List(%))
+
+    s1 + s2 ==
+      s1 = 0 => s2
+      s2 = 0 => s1
+      conv [convert("+"::Symbol), s1, s2]$List(%)
+
+    s1 * s2 ==
+      s1 = 0 or s2 = 0 => 0
+      s1 = 1 => s2
+      s2 = 1 => s1
+      conv [convert("*"::Symbol), s1, s2]$List(%)
+
+    s1:% ** n:Integer ==
+      s1 = 0 and n > 0 => 0
+      s1 = 1 or zero? n => 1
+--      one? n => s1
+      (n = 1) => s1
+      conv [convert("**"::Symbol), s1, convert n]$List(%)
+
+    s1:% ** n:NonNegativeInteger == s1 ** (n::Integer)
+
+    s1 / s2 ==
+      s2 = 1 => s1
+      conv [convert("/"::Symbol), s1, s2]$List(%)
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain INT Integer}
 The function {\bf one?} has been rewritten back to its original form.
 The NAG version called a lisp primitive that exists only in Codemist
@@ -32279,24 +32774,6 @@ leq
 --E 16
 )spool
  
-Dx: LODO(EXPR INT, f +-> D(f, x))
-Dx := D()
-Dop:= Dx^3 + G/x^2*Dx + H/x^3 - 1
-n == 3
-phi == reduce(+,[subscript(s,[i])*exp(x)/x^i for i in 0..n])
-phi1 ==  Dop(phi) / exp x
-phi2 == phi1 *x**(n+3)
-phi3 == retract(phi2)@(POLY INT)
-pans == phi3 ::UP(x,POLY INT)
-pans1 == [coefficient(pans, (n+3-i) :: NNI) for i in 2..n+1]
-leq == solve(pans1,[subscript(s,[i]) for i in 1..n])
-leq
-n==4
-leq
-n==7
-leq
-)spool
-)lisp (bye)
 @
 <<LinearOrdinaryDifferentialOperator.help>>=
 ====================================================================
@@ -32499,6 +32976,9 @@ o $AXIOM/doc/src/algebra/lodo.spad.dvi
 @
 \pagehead{LinearOrdinaryDifferentialOperator}{LODO}
 \pagepic{ps/v103linearordinarydifferentialoperator.ps}{LODO}{1.00}
+See also:\\
+\refto{LinearOrdinaryDifferentialOperator1}{LODO1}
+\refto{LinearOrdinaryDifferentialOperator2}{LODO2}
 <<domain LODO LinearOrdinaryDifferentialOperator>>=
 )abbrev domain LODO LinearOrdinaryDifferentialOperator
 ++ Author: Manuel Bronstein
@@ -32925,6 +33405,9 @@ o $AXIOM/doc/src/algebra/lodo.spad.dvi
 @
 \pagehead{LinearOrdinaryDifferentialOperator1}{LODO1}
 \pagepic{ps/v103linearordinarydifferentialoperator1.ps}{LODO1}{1.00}
+See also:\\
+\refto{LinearOrdinaryDifferentialOperator}{LODO}
+\refto{LinearOrdinaryDifferentialOperator2}{LODO2}
 <<domain LODO1 LinearOrdinaryDifferentialOperator1>>=
 )abbrev domain LODO1 LinearOrdinaryDifferentialOperator1
 ++ Author: Manuel Bronstein
@@ -33465,6 +33948,9 @@ o $AXIOM/doc/src/algebra/lodo.spad.dvi
 @
 \pagehead{LinearOrdinaryDifferentialOperator2}{LODO2}
 \pagepic{ps/v103linearordinarydifferentialoperator2.ps}{LODO2}{1.00}
+See also:\\
+\refto{LinearOrdinaryDifferentialOperator}{LODO}
+\refto{LinearOrdinaryDifferentialOperator1}{LODO1}
 <<domain LODO2 LinearOrdinaryDifferentialOperator2>>=
 )abbrev domain LODO2 LinearOrdinaryDifferentialOperator2
 ++ Author: Stephen M. Watt, Manuel Bronstein
@@ -35085,6 +35571,2600 @@ MakeCachableSet(S:SetCategory): Exports == 
Implementation where
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MATRIX Matrix}
+<<Matrix.input>>=
+-- matrix.spad.pamphlet Matrix.input
+)spool Matrix.output
+)set message test on
+)set message auto off
+)clear all
+--S 1 of 38
+m : Matrix(Integer) := new(3,3,0)
+--R 
+--R
+--R        +0  0  0+
+--R        |       |
+--R   (1)  |0  0  0|
+--R        |       |
+--R        +0  0  0+
+--R                                                         Type: Matrix 
Integer
+--E 1
+
+--S 2 of 38
+setelt(m,2,3,5)
+--R 
+--R
+--R   (2)  5
+--R                                                        Type: 
PositiveInteger
+--E 2
+
+--S 3 of 38
+m(1,2) := 10
+--R 
+--R
+--R   (3)  10
+--R                                                        Type: 
PositiveInteger
+--E 3
+
+--S 4 of 38
+m
+--R 
+--R
+--R        +0  10  0+
+--R        |        |
+--R   (4)  |0  0   5|
+--R        |        |
+--R        +0  0   0+
+--R                                                         Type: Matrix 
Integer
+--E 4
+
+--S 5 of 38
+matrix [ [1,2,3,4],[0,9,8,7] ]
+--R 
+--R
+--R        +1  2  3  4+
+--R   (5)  |          |
+--R        +0  9  8  7+
+--R                                                         Type: Matrix 
Integer
+--E 5
+
+--S 6 of 38
+dm := diagonalMatrix [1,x**2,x**3,x**4,x**5]
+--R 
+--R
+--R        +1  0   0   0   0 +
+--R        |                 |
+--R        |    2            |
+--R        |0  x   0   0   0 |
+--R        |                 |
+--R        |        3        |
+--R   (6)  |0  0   x   0   0 |
+--R        |                 |
+--R        |            4    |
+--R        |0  0   0   x   0 |
+--R        |                 |
+--R        |                5|
+--R        +0  0   0   0   x +
+--R                                              Type: Matrix Polynomial 
Integer
+--E 6
+
+--S 7 of 38
+setRow!(dm,5,vector [1,1,1,1,1])
+--R 
+--R
+--R        +1  0   0   0   0+
+--R        |                |
+--R        |    2           |
+--R        |0  x   0   0   0|
+--R        |                |
+--R   (7)  |        3       |
+--R        |0  0   x   0   0|
+--R        |                |
+--R        |            4   |
+--R        |0  0   0   x   0|
+--R        |                |
+--R        +1  1   1   1   1+
+--R                                              Type: Matrix Polynomial 
Integer
+--E 7
+
+--S 8 of 38
+setColumn!(dm,2,vector [y,y,y,y,y])
+--R 
+--R
+--R        +1  y  0   0   0+
+--R        |               |
+--R        |0  y  0   0   0|
+--R        |               |
+--R        |       3       |
+--R   (8)  |0  y  x   0   0|
+--R        |               |
+--R        |           4   |
+--R        |0  y  0   x   0|
+--R        |               |
+--R        +1  y  1   1   1+
+--R                                              Type: Matrix Polynomial 
Integer
+--E 8
+
+--S 9 of 38
+cdm := copy(dm)
+--R 
+--R
+--R        +1  y  0   0   0+
+--R        |               |
+--R        |0  y  0   0   0|
+--R        |               |
+--R        |       3       |
+--R   (9)  |0  y  x   0   0|
+--R        |               |
+--R        |           4   |
+--R        |0  y  0   x   0|
+--R        |               |
+--R        +1  y  1   1   1+
+--R                                              Type: Matrix Polynomial 
Integer
+--E 9
+
+--S 10 of 38
+setelt(dm,4,1,1-x**7)
+--R 
+--R
+--R            7
+--R   (10)  - x  + 1
+--R                                                     Type: Polynomial 
Integer
+--E 10
+
+--S 11 of 38
+[dm,cdm]
+--R 
+--R
+--R          +   1      y  0   0   0+ +1  y  0   0   0+
+--R          |                      | |               |
+--R          |   0      y  0   0   0| |0  y  0   0   0|
+--R          |                      | |               |
+--R          |              3       | |       3       |
+--R   (11)  [|   0      y  x   0   0|,|0  y  x   0   0|]
+--R          |                      | |               |
+--R          |   7              4   | |           4   |
+--R          |- x  + 1  y  0   x   0| |0  y  0   x   0|
+--R          |                      | |               |
+--R          +   1      y  1   1   1+ +1  y  1   1   1+
+--R                                         Type: List Matrix Polynomial 
Integer
+--E 11
+
+--S 12 of 38
+subMatrix(dm,2,3,2,4)
+--R 
+--R
+--R         +y  0   0+
+--R   (12)  |        |
+--R         |    3   |
+--R         +y  x   0+
+--R                                              Type: Matrix Polynomial 
Integer
+--E 12
+
+--S 13 of 38
+d := diagonalMatrix [1.2,-1.3,1.4,-1.5]
+--R 
+--R
+--R         +1.2   0.0   0.0   0.0 +
+--R         |                      |
+--R         |0.0  - 1.3  0.0   0.0 |
+--R   (13)  |                      |
+--R         |0.0   0.0   1.4   0.0 |
+--R         |                      |
+--R         +0.0   0.0   0.0  - 1.5+
+--R                                                           Type: Matrix 
Float
+--E 13
+
+--S 14 of 38
+e := matrix [ [6.7,9.11],[-31.33,67.19] ]
+--R 
+--R
+--R         +  6.7    9.11 +
+--R   (14)  |              |
+--R         +- 31.33  67.19+
+--R                                                           Type: Matrix 
Float
+--E 14
+
+--S 15 of 38
+setsubMatrix!(d,1,2,e)
+--R 
+--R
+--R         +1.2    6.7    9.11    0.0 +
+--R         |                          |
+--R         |0.0  - 31.33  67.19   0.0 |
+--R   (15)  |                          |
+--R         |0.0    0.0     1.4    0.0 |
+--R         |                          |
+--R         +0.0    0.0     0.0   - 1.5+
+--R                                                           Type: Matrix 
Float
+--E 15
+
+--S 16 of 38
+d
+--R 
+--R
+--R         +1.2    6.7    9.11    0.0 +
+--R         |                          |
+--R         |0.0  - 31.33  67.19   0.0 |
+--R   (16)  |                          |
+--R         |0.0    0.0     1.4    0.0 |
+--R         |                          |
+--R         +0.0    0.0     0.0   - 1.5+
+--R                                                           Type: Matrix 
Float
+--E 16
+
+--S 17 of 38
+a := matrix [ [1/2,1/3,1/4],[1/5,1/6,1/7] ]
+--R 
+--R
+--R         +1  1  1+
+--R         |-  -  -|
+--R         |2  3  4|
+--R   (17)  |       |
+--R         |1  1  1|
+--R         |-  -  -|
+--R         +5  6  7+
+--R                                                Type: Matrix Fraction 
Integer
+--E 17
+
+--S 18 of 38
+b := matrix [ [3/5,3/7,3/11],[3/13,3/17,3/19] ] 
+--R 
+--R
+--R         +3   3    3+
+--R         |-   -   --|
+--R         |5   7   11|
+--R   (18)  |          |
+--R         | 3   3   3|
+--R         |--  --  --|
+--R         +13  17  19+
+--R                                                Type: Matrix Fraction 
Integer
+--E 18
+
+--S 19 of 38
+horizConcat(a,b)
+--R 
+--R
+--R         +1  1  1  3   3    3+
+--R         |-  -  -  -   -   --|
+--R         |2  3  4  5   7   11|
+--R   (19)  |                   |
+--R         |1  1  1   3   3   3|
+--R         |-  -  -  --  --  --|
+--R         +5  6  7  13  17  19+
+--R                                                Type: Matrix Fraction 
Integer
+--E 19
+
+--S 20 of 38
+vab := vertConcat(a,b)
+--R 
+--R
+--R         +1   1   1 +
+--R         |-   -   - |
+--R         |2   3   4 |
+--R         |          |
+--R         |1   1   1 |
+--R         |-   -   - |
+--R         |5   6   7 |
+--R   (20)  |          |
+--R         |3   3    3|
+--R         |-   -   --|
+--R         |5   7   11|
+--R         |          |
+--R         | 3   3   3|
+--R         |--  --  --|
+--R         +13  17  19+
+--R                                                Type: Matrix Fraction 
Integer
+--E 20
+
+--S 21 of 38
+transpose vab
+--R 
+--R
+--R         +1  1  3    3+
+--R         |-  -  -   --|
+--R         |2  5  5   13|
+--R         |            |
+--R         |1  1  3    3|
+--R   (21)  |-  -  -   --|
+--R         |3  6  7   17|
+--R         |            |
+--R         |1  1   3   3|
+--R         |-  -  --  --|
+--R         +4  7  11  19+
+--R                                                Type: Matrix Fraction 
Integer
+--E 21
+
+--S 22 of 38
+m := matrix [ [1,2],[3,4] ]
+--R 
+--R
+--R         +1  2+
+--R   (22)  |    |
+--R         +3  4+
+--R                                                         Type: Matrix 
Integer
+--E 22
+
+--S 23 of 38
+4 * m * (-5)
+--R 
+--R
+--R         +- 20  - 40+
+--R   (23)  |          |
+--R         +- 60  - 80+
+--R                                                         Type: Matrix 
Integer
+--E 23
+
+--S 24 of 38
+n := matrix([ [1,0,-2],[-3,5,1] ])
+--R 
+--R
+--R         + 1   0  - 2+
+--R   (24)  |           |
+--R         +- 3  5   1 +
+--R                                                         Type: Matrix 
Integer
+--E 24
+
+--S 25 of 38
+m * n
+--R 
+--R
+--R         +- 5  10   0 +
+--R   (25)  |            |
+--R         +- 9  20  - 2+
+--R                                                         Type: Matrix 
Integer
+--E 25
+
+--S 26 of 38
+vec := column(n,3)
+--R 
+--R
+--R   (26)  [- 2,1]
+--R                                                         Type: Vector 
Integer
+--E 26
+
+--S 27 of 38
+vec * m
+--R 
+--R
+--R   (27)  [1,0]
+--R                                                         Type: Vector 
Integer
+--E 27
+
+--S 28 of 38
+m * vec
+--R 
+--R
+--R   (28)  [0,- 2]
+--R                                                         Type: Vector 
Integer
+--E 28
+
+--S 29 of 38
+hilb := matrix([ [1/(i + j) for i in 1..3] for j in 1..3])
+--R 
+--R
+--R         +1  1  1+
+--R         |-  -  -|
+--R         |2  3  4|
+--R         |       |
+--R         |1  1  1|
+--R   (29)  |-  -  -|
+--R         |3  4  5|
+--R         |       |
+--R         |1  1  1|
+--R         |-  -  -|
+--R         +4  5  6+
+--R                                                Type: Matrix Fraction 
Integer
+--E 29
+
+--S 30 of 38
+inverse(hilb)
+--R 
+--R
+--R         + 72    - 240   180 +
+--R         |                   |
+--R   (30)  |- 240   900   - 720|
+--R         |                   |
+--R         + 180   - 720   600 +
+--R                                     Type: Union(Matrix Fraction 
Integer,...)
+--E 30
+
+--S 31 of 38
+mm := matrix([ [1,2,3,4], [5,6,7,8], [9,10,11,12], [13,14,15,16] ])
+--R 
+--R
+--R         +1   2   3   4 +
+--R         |              |
+--R         |5   6   7   8 |
+--R   (31)  |              |
+--R         |9   10  11  12|
+--R         |              |
+--R         +13  14  15  16+
+--R                                                         Type: Matrix 
Integer
+--E 31
+
+--S 32 of 38
+inverse(mm)
+--R 
+--R
+--R   (32)  "failed"
+--R                                                    Type: 
Union("failed",...)
+--E 32
+
+--S 33 of 38
+determinant(mm)
+--R 
+--R
+--R   (33)  0
+--R                                                     Type: 
NonNegativeInteger
+--E 33
+
+--S 34 of 38
+trace(mm)
+--R 
+--R
+--R   (34)  34
+--R                                                        Type: 
PositiveInteger
+--E 34
+
+--S 35 of 38
+rank(mm)
+--R 
+--R
+--R   (35)  2
+--R                                                        Type: 
PositiveInteger
+--E 35
+
+--S 36 of 38
+nullity(mm)
+--R 
+--R
+--R   (36)  2
+--R                                                        Type: 
PositiveInteger
+--E 36
+
+--S 37 of 38
+nullSpace(mm)
+--R 
+--R
+--R   (37)  [[1,- 2,1,0],[2,- 3,0,1]]
+--R                                                    Type: List Vector 
Integer
+--E 37
+
+--S 38 of 38
+rowEchelon(mm)
+--R 
+--R
+--R         +1  2  3  4 +
+--R         |           |
+--R         |0  4  8  12|
+--R   (38)  |           |
+--R         |0  0  0  0 |
+--R         |           |
+--R         +0  0  0  0 +
+--R                                                         Type: Matrix 
Integer
+--E 38
+)spool
+)lisp (bye)
+@
+<<Matrix.help>>=
+====================================================================
+Matrix examples
+====================================================================
+
+The Matrix domain provides arithmetic operations on matrices
+and standard functions from linear algebra.
+This domain is similar to the TwoDimensionalArray domain, except
+that the entries for Matrix must belong to a  Ring.
+
+====================================================================
+Creating Matrices
+====================================================================
+
+There are many ways to create a matrix from a collection of values or
+from existing matrices.
+
+If the matrix has almost all items equal to the same value, use new to
+create a matrix filled with that value and then reset the entries that
+are different.
+
+  m : Matrix(Integer) := new(3,3,0)
+    +0  0  0+
+    |       |
+    |0  0  0|
+    |       |
+    +0  0  0+
+                      Type: Matrix Integer
+
+To change the entry in the second row, third column to 5, use setelt.
+
+  setelt(m,2,3,5)
+    5
+                      Type: PositiveInteger
+
+An alternative syntax is to use assignment.
+
+  m(1,2) := 10
+    10
+                      Type: PositiveInteger
+
+The matrix was destructively modified.
+
+  m
+    +0  10  0+
+    |        |
+    |0  0   5|
+    |        |
+    +0  0   0+
+                      Type: Matrix Integer
+
+If you already have the matrix entries as a list of lists, use matrix.
+
+  matrix [ [1,2,3,4],[0,9,8,7] ]
+    +1  2  3  4+
+    |          |
+    +0  9  8  7+
+                      Type: Matrix Integer
+
+If the matrix is diagonal, use diagonalMatrix.
+
+  dm := diagonalMatrix [1,x**2,x**3,x**4,x**5]
+        +1  0   0   0   0 +
+        |                 |
+        |    2            |
+        |0  x   0   0   0 |
+        |                 |
+        |        3        |
+        |0  0   x   0   0 |
+        |                 |
+        |            4    |
+        |0  0   0   x   0 |
+        |                 |
+        |                5|
+        +0  0   0   0   x +
+                     Type: Matrix Polynomial Integer
+
+Use setRow and setColumn to change a row or column of a matrix.
+
+  setRow!(dm,5,vector [1,1,1,1,1])
+        +1  0   0   0   0+
+        |                |
+        |    2           |
+        |0  x   0   0   0|
+        |                |
+        |        3       |
+        |0  0   x   0   0|
+        |                |
+        |            4   |
+        |0  0   0   x   0|
+        |                |
+        +1  1   1   1   1+
+                    Type: Matrix Polynomial Integer
+
+  setColumn!(dm,2,vector [y,y,y,y,y])
+        +1  y  0   0   0+
+        |               |
+        |0  y  0   0   0|
+        |               |
+        |       3       |
+        |0  y  x   0   0|
+        |               |
+        |           4   |
+        |0  y  0   x   0|
+        |               |
+        +1  y  1   1   1+
+                    Type: Matrix Polynomial Integer
+
+Use copy to make a copy of a matrix.
+
+  cdm := copy(dm)
+        +1  y  0   0   0+
+        |               |
+        |0  y  0   0   0|
+        |               |
+        |       3       |
+        |0  y  x   0   0|
+        |               |
+        |           4   |
+        |0  y  0   x   0|
+        |               |
+        +1  y  1   1   1+
+                    Type: Matrix Polynomial Integer
+
+This is useful if you intend to modify a matrix destructively but
+want a copy of the original.
+
+  setelt(dm,4,1,1-x**7)
+        7
+     - x  + 1
+                    Type: Polynomial Integer
+
+  [dm,cdm]
+          +   1      y  0   0   0+ +1  y  0   0   0+
+          |                      | |               |
+          |   0      y  0   0   0| |0  y  0   0   0|
+          |                      | |               |
+          |              3       | |       3       |
+         [|   0      y  x   0   0|,|0  y  x   0   0|]
+          |                      | |               |
+          |   7              4   | |           4   |
+          |- x  + 1  y  0   x   0| |0  y  0   x   0|
+          |                      | |               |
+          +   1      y  1   1   1+ +1  y  1   1   1+
+                     Type: List Matrix Polynomial Integer
+
+Use subMatrix to extract part of an existing matrix.  The syntax is 
+subMatrix(m, firstrow, lastrow, firstcol, lastcol).
+
+  subMatrix(dm,2,3,2,4)
+         +y  0   0+
+         |        |
+         |    3   |
+         +y  x   0+
+                     Type: Matrix Polynomial Integer
+
+To change a submatrix, use setsubMatrix.
+
+  d := diagonalMatrix [1.2,-1.3,1.4,-1.5]
+         +1.2   0.0   0.0   0.0 +
+         |                      |
+         |0.0  - 1.3  0.0   0.0 |
+         |                      |
+         |0.0   0.0   1.4   0.0 |
+         |                      |
+         +0.0   0.0   0.0  - 1.5+
+                     Type: Matrix Float
+
+If e is too big to fit where you specify, an error message is
+displayed.  Use subMatrix to extract part of e, if necessary.
+
+  e := matrix [ [6.7,9.11],[-31.33,67.19] ]
+         +  6.7    9.11 +
+         |              |
+         +- 31.33  67.19+
+                      Type: Matrix Float
+
+This changes the submatrix of d whose upper left corner is at the
+first row and second column and whose size is that of e.
+
+  setsubMatrix!(d,1,2,e)
+         +1.2    6.7    9.11    0.0 +
+         |                          |
+         |0.0  - 31.33  67.19   0.0 |
+         |                          |
+         |0.0    0.0     1.4    0.0 |
+         |                          |
+         +0.0    0.0     0.0   - 1.5+
+                       Type: Matrix Float
+
+  d
+         +1.2    6.7    9.11    0.0 +
+         |                          |
+         |0.0  - 31.33  67.19   0.0 |
+         |                          |
+         |0.0    0.0     1.4    0.0 |
+         |                          |
+         +0.0    0.0     0.0   - 1.5+
+                        Type: Matrix Float
+
+Matrices can be joined either horizontally or vertically to make
+new matrices.
+
+  a := matrix [ [1/2,1/3,1/4],[1/5,1/6,1/7] ]
+         +1  1  1+
+         |-  -  -|
+         |2  3  4|
+         |       |
+         |1  1  1|
+         |-  -  -|
+         +5  6  7+
+                         Type: Matrix Fraction Integer
+
+  b := matrix [ [3/5,3/7,3/11],[3/13,3/17,3/19] ] 
+         +3   3    3+
+         |-   -   --|
+         |5   7   11|
+         |          |
+         | 3   3   3|
+         |--  --  --|
+         +13  17  19+
+                         Type: Matrix Fraction Integer
+
+Use horizConcat to append them side to side.  The two matrices must
+have the same number of rows.
+
+  horizConcat(a,b)
+         +1  1  1  3   3    3+
+         |-  -  -  -   -   --|
+         |2  3  4  5   7   11|
+         |                   |
+         |1  1  1   3   3   3|
+         |-  -  -  --  --  --|
+         +5  6  7  13  17  19+
+                         Type: Matrix Fraction Integer
+
+Use vertConcat to stack one upon the other.  The two matrices must
+have the same number of columns.
+
+  vab := vertConcat(a,b)
+         +1   1   1 +
+         |-   -   - |
+         |2   3   4 |
+         |          |
+         |1   1   1 |
+         |-   -   - |
+         |5   6   7 |
+         |          |
+         |3   3    3|
+         |-   -   --|
+         |5   7   11|
+         |          |
+         | 3   3   3|
+         |--  --  --|
+         +13  17  19+
+                         Type: Matrix Fraction Integer
+
+The operation transpose is used to create a new matrix by reflection
+across the main diagonal.
+
+  transpose vab
+         +1  1  3    3+
+         |-  -  -   --|
+         |2  5  5   13|
+         |            |
+         |1  1  3    3|
+         |-  -  -   --|
+         |3  6  7   17|
+         |            |
+         |1  1   3   3|
+         |-  -  --  --|
+         +4  7  11  19+
+                         Type: Matrix Fraction Integer
+
+====================================================================
+Operations on Matrices
+====================================================================
+
+Axiom provides both left and right scalar multiplication.
+
+  m := matrix [ [1,2],[3,4] ]
+         +1  2+
+         |    |
+         +3  4+
+                          Type: Matrix Integer
+
+  4 * m * (-5)
+         +- 20  - 40+
+         |          |
+         +- 60  - 80+
+                          Type: Matrix Integer
+
+You can add, subtract, and multiply matrices provided, of course, that
+the matrices have compatible dimensions.  If not, an error message is
+displayed.
+
+  n := matrix([ [1,0,-2],[-3,5,1] ])
+         + 1   0  - 2+
+         |           |
+         +- 3  5   1 +
+                          Type: Matrix Integer
+
+This following product is defined but n * m is not.
+
+  m * n
+         +- 5  10   0 +
+         |            |
+         +- 9  20  - 2+
+                          Type: Matrix Integer
+
+The operations nrows and ncols return the number of rows and columns
+of a matrix.  You can extract a row or a column of a matrix using the
+operations row and column.  The object returned is a Vector.
+
+Here is the third column of the matrix n.
+
+  vec := column(n,3)
+     [- 2,1]
+                          Type: Vector Integer
+
+You can multiply a matrix on the left by a "row vector" and on the right
+by a "column vector".
+
+  vec * m
+     [1,0]
+                          Type: Vector Integer
+
+Of course, the dimensions of the vector and the matrix must be compatible
+or an error message is returned.
+
+  m * vec
+    [0,- 2]
+                          Type: Vector Integer
+
+The operation inverse computes the inverse of a matrix if the matrix
+is invertible, and returns "failed" if not.
+
+This Hilbert matrix is invertible.
+
+  hilb := matrix([ [1/(i + j) for i in 1..3] for j in 1..3])
+         +1  1  1+
+         |-  -  -|
+         |2  3  4|
+         |       |
+         |1  1  1|
+         |-  -  -|
+         |3  4  5|
+         |       |
+         |1  1  1|
+         |-  -  -|
+         +4  5  6+
+                          Type: Matrix Fraction Integer
+
+  inverse(hilb)
+         + 72    - 240   180 +
+         |                   |
+         |- 240   900   - 720|
+         |                   |
+         + 180   - 720   600 +
+                          Type: Union(Matrix Fraction Integer,...)
+
+This matrix is not invertible.
+
+  mm := matrix([ [1,2,3,4], [5,6,7,8], [9,10,11,12], [13,14,15,16] ])
+         +1   2   3   4 +
+         |              |
+         |5   6   7   8 |
+         |              |
+         |9   10  11  12|
+         |              |
+         +13  14  15  16+
+                           Type: Matrix Integer
+
+  inverse(mm)
+     "failed"
+                           Type: Union("failed",...)
+
+The operation determinant computes the determinant of a matrix
+provided that the entries of the matrix belong to a CommutativeRing.
+
+The above matrix mm is not invertible and, hence, must have determinant 0.
+
+  determinant(mm)
+    0
+                           Type: NonNegativeInteger
+
+The operation trace computes the trace of a square matrix.
+
+  trace(mm)
+    34
+                           Type: PositiveInteger
+
+The operation rank computes the rank of a matrix: the maximal number
+of linearly independent rows or columns.
+
+  rank(mm)
+    2
+                           Type: PositiveInteger
+
+The operation nullity computes the nullity of a matrix: the dimension
+of its null space.
+
+  nullity(mm)
+    2
+                           Type: PositiveInteger
+
+The operation nullSpace returns a list containing a basis for the null
+space of a matrix.  Note that the nullity is the number of elements in
+a basis for the null space.
+
+  nullSpace(mm)
+    [[1,- 2,1,0],[2,- 3,0,1]]
+                           Type: List Vector Integer
+
+The operation rowEchelon returns the row echelon form of a matrix.  It
+is easy to see that the rank of this matrix is two and that its
+nullity is also two.
+
+  rowEchelon(mm)
+         +1  2  3  4 +
+         |           |
+         |0  4  8  12|
+         |           |
+         |0  0  0  0 |
+         |           |
+         +0  0  0  0 +
+                           Type: Matrix Integer
+
+See Also
+o )help OneDimensionalArray
+o )help TwoDimensionalArray
+o )help Vector
+o )help Permanent
+o )show Matrix
+o $AXIOM/doc/src/algebra/matrix.spad.dvi
+
+@
+\pagehead{Matrix}{MATRIX}
+\pagepic{ps/v103matrix.ps}{MATRIX}{1.00}
+See also:\\
+\refto{IndexedMatrix}{IMATRIX}
+\refto{RectangularMatrix}{RMATRIX}
+\refto{SquareMatrix}{SQMATRIX}
+<<domain MATRIX Matrix>>=
+)abbrev domain MATRIX Matrix
+++ Author: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: IndexedMatrix, RectangularMatrix, SquareMatrix
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{Matrix} is a matrix domain where 1-based indexing is used
+++   for both rows and columns.
+Matrix(R): Exports == Implementation where
+  R : Ring
+  Row ==> Vector R
+  Col ==> Vector R
+  mnRow ==> 1
+  mnCol ==> 1
+  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,$)
+  MATSTOR ==> StorageEfficientMatrixOperations(R)
+ 
+  Exports ==> MatrixCategory(R,Row,Col) with
+    diagonalMatrix: Vector R -> $
+      ++ \spad{diagonalMatrix(v)} returns a diagonal matrix where the elements
+      ++ of v appear on the diagonal.
+
+    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
+
+    if R has Field then
+      inverse: $ -> Union($,"failed")
+        ++ \spad{inverse(m)} returns the inverse of the matrix m. 
+        ++ If the matrix is not invertible, "failed" is returned.
+        ++ Error: if the matrix is not square.
+--     matrix: Vector Vector R -> $
+--       ++ \spad{matrix(v)} converts the vector of vectors v to a matrix, 
where
+--       ++ the vector of vectors is viewed as a vector of the rows of the
+--       ++ matrix
+--     diagonalMatrix: Vector $ -> $
+--       ++ \spad{diagonalMatrix([m1,...,mk])} creates a block diagonal matrix
+--       ++ M with block matrices {\em m1},...,{\em mk} down the diagonal,
+--       ++ with 0 block matrices elsewhere.
+--     vectorOfVectors: $ -> Vector Vector R
+--       ++ \spad{vectorOfVectors(m)} returns the rows of the matrix m as a
+--       ++ vector of vectors
+ 
+  Implementation ==>
+   InnerIndexedTwoDimensionalArray(R,mnRow,mnCol,Row,Col) add
+    minr ==> minRowIndex
+    maxr ==> maxRowIndex
+    minc ==> minColIndex
+    maxc ==> maxColIndex
+    mini ==> minIndex
+    maxi ==> maxIndex
+ 
+    minRowIndex x == mnRow
+    minColIndex x == mnCol
+ 
+    swapRows_!(x,i1,i2) ==
+        (i1 < minRowIndex(x)) or (i1 > maxRowIndex(x)) or _
+           (i2 < minRowIndex(x)) or (i2 > maxRowIndex(x)) =>
+             error "swapRows!: index out of range"
+        i1 = i2 => x
+        minRow := minRowIndex x
+        xx := x pretend PrimitiveArray(PrimitiveArray(R))
+        n1 := i1 - minRow; n2 := i2 - minRow
+        row1 := qelt(xx,n1)
+        qsetelt_!(xx,n1,qelt(xx,n2))
+        qsetelt_!(xx,n2,row1)
+        xx pretend $
+ 
+    positivePower:($,Integer,NonNegativeInteger) -> $
+    positivePower(x,n,nn) ==
+--      one? n => x
+      (n = 1) => x
+      -- no need to allocate space for 3 additional matrices
+      n = 2 => x * x
+      n = 3 => x * x * x
+      n = 4 => (y := x * x; y * y)
+      a := new(nn,nn,0) pretend Matrix(R)
+      b := new(nn,nn,0) pretend Matrix(R)
+      c := new(nn,nn,0) pretend Matrix(R)
+      xx := x pretend Matrix(R)
+      power_!(a,b,c,xx,n :: NonNegativeInteger)$MATSTOR pretend $
+ 
+    x:$ ** n:NonNegativeInteger ==
+      not((nn := nrows x) = ncols x) =>
+        error "**: matrix must be square"
+      zero? n => scalarMatrix(nn,1)
+      positivePower(x,n,nn)
+ 
+    if R has commutative("*") then
+ 
+        determinant x == determinant(x)$MATLIN
+        minordet    x == minordet(x)$MATLIN
+ 
+    if R has EuclideanDomain then
+ 
+        rowEchelon  x == rowEchelon(x)$MATLIN
+ 
+    if R has IntegralDomain then
+ 
+        rank        x == rank(x)$MATLIN
+        nullity     x == nullity(x)$MATLIN
+        nullSpace   x == nullSpace(x)$MATLIN
+ 
+    if R has Field then
+ 
+        inverse     x == inverse(x)$MATLIN
+ 
+        x:$ ** n:Integer ==
+          nn := nrows x
+          not(nn = ncols x) =>
+            error "**: matrix must be square"
+          zero? n => scalarMatrix(nn,1)
+          positive? n => positivePower(x,n,nn)
+          (xInv := inverse x) case "failed" =>
+            error "**: matrix must be invertible"
+          positivePower(xInv :: $,-n,nn)
+ 
+--     matrix(v: Vector Vector R) ==
+--       (rows := # v) = 0 => new(0,0,0)
+--       -- error check: this is a top level function
+--       cols := # v.mini(v)
+--       for k in (mini(v) + 1)..maxi(v) repeat
+--         cols ^= # v.k => error "matrix: rows of different lengths"
+--       ans := new(rows,cols,0)
+--       for i in minr(ans)..maxr(ans) for k in mini(v)..maxi(v) repeat
+--         vv := v.k
+--         for j in minc(ans)..maxc(ans) for l in mini(vv)..maxi(vv) repeat
+--           ans(i,j) := vv.l
+--       ans
+ 
+    diagonalMatrix(v: Vector R) ==
+      n := #v; ans := zero(n,n)
+      for i in minr(ans)..maxr(ans) for j in minc(ans)..maxc(ans) _
+          for k in mini(v)..maxi(v) repeat qsetelt_!(ans,i,j,qelt(v,k))
+      ans
+ 
+--     diagonalMatrix(vec: Vector $) ==
+--       rows : NonNegativeInteger := 0
+--       cols : NonNegativeInteger := 0
+--       for r in mini(vec)..maxi(vec) repeat
+--         mat := vec.r
+--         rows := rows + nrows mat; cols := cols + ncols mat
+--       ans := zero(rows,cols)
+--       loR := minr ans; loC := minc ans
+--       for r in mini(vec)..maxi(vec) repeat
+--         mat := vec.r
+--         hiR := loR + nrows(mat) - 1; hiC := loC + nrows(mat) - 1
+--         for i in loR..hiR for k in minr(mat)..maxr(mat) repeat
+--           for j in loC..hiC for l in minc(mat)..maxc(mat) repeat
+--             ans(i,j) := mat(k,l)
+--         loR := hiR + 1; loC := hiC + 1
+--       ans
+ 
+--     vectorOfVectors x ==
+--       vv : Vector Vector R := new(nrows x,0)
+--       cols := ncols x
+--       for k in mini(vv)..maxi(vv) repeat
+--         vv.k := new(cols,0)
+--       for i in minr(x)..maxr(x) for k in mini(vv)..maxi(vv) repeat
+--         v := vv.k
+--         for j in minc(x)..maxc(x) for l in mini(v)..maxi(v) repeat
+--           v.l := x(i,j)
+--       vv
+ 
+    if R has ConvertibleTo InputForm then
+      convert(x:$):InputForm ==
+         convert [convert("matrix"::Symbol)@InputForm,
+                  convert listOfLists x]$List(InputForm)
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MODMON ModMonic}
+\pagehead{ModMonic}{MODMON}
+\pagepic{ps/v103modmonic.ps}{MODMON}{1.00}
+<<domain MODMON ModMonic>>=
+)abbrev domain MODMON ModMonic
+++ Description:
+++ This package \undocumented
+-- following line prevents caching ModMonic
+)bo PUSH('ModMonic, $mutableDomains)
+ 
+ModMonic(R,Rep): C == T
+ where
+  R: Ring
+  Rep: UnivariatePolynomialCategory(R)
+  C == UnivariatePolynomialCategory(R) with
+  --operations
+    setPoly : Rep -> Rep
+       ++ setPoly(x) \undocumented
+    modulus : -> Rep
+       ++ modulus() \undocumented
+    reduce: Rep -> %
+       ++ reduce(x) \undocumented
+    lift: % -> Rep --reduce lift = identity
+       ++ lift(x) \undocumented
+    coerce: Rep -> %
+       ++ coerce(x) \undocumented
+    Vectorise: % -> Vector(R)
+       ++ Vectorise(x) \undocumented
+    UnVectorise: Vector(R) -> %
+       ++ UnVectorise(v) \undocumented
+    An: % -> Vector(R)
+       ++ An(x) \undocumented
+    pow : -> PrimitiveArray(%)
+       ++ pow() \undocumented
+    computePowers : -> PrimitiveArray(%)
+       ++ computePowers() \undocumented
+    if R has FiniteFieldCategory then
+       frobenius: % -> %
+       ++ frobenius(x) \undocumented
+    --LinearTransf: (%,Vector(R)) -> SquareMatrix<deg> R
+  --assertions
+    if R has Finite then Finite
+  T == add
+    --constants
+      m:Rep := monomial(1,1)$Rep --| degree(m) > 0 and LeadingCoef(m) = R$1
+      d := degree(m)$Rep
+      d1 := (d-1):NonNegativeInteger
+      twod := 2*d1+1
+      frobenius?:Boolean := R has FiniteFieldCategory
+      --VectorRep:= DirectProduct(d:NonNegativeInteger,R)
+    --declarations
+      x,y: %
+      p: Rep
+      d,n: Integer
+      e,k1,k2: NonNegativeInteger
+      c: R
+      --vect: Vector(R)
+      power:PrimitiveArray(%)
+      frobeniusPower:PrimitiveArray(%)
+      computeFrobeniusPowers : () -> PrimitiveArray(%)
+    --representations
+    --mutable m    --take this out??
+    --define
+      power := new(0,0)
+      frobeniusPower := new(0,0)
+      setPoly (mon : Rep) ==
+        mon =$Rep m => mon
+        oldm := m
+        leadingCoefficient mon ^= 1 => error "polynomial must be monic"
+        -- following copy code needed since FFPOLY can modify mon
+        copymon:Rep:= 0
+        while not zero? mon repeat
+           copymon := monomial(leadingCoefficient mon, degree mon)$Rep + 
copymon
+           mon := reductum mon
+        m := copymon
+        d := degree(m)$Rep
+        d1 := (d-1)::NonNegativeInteger
+        twod := 2*d1+1
+        power := computePowers()
+        if frobenius? then
+          degree(oldm)>1 and not((oldm exquo$Rep m) case "failed") =>
+              for i in 1..d1 repeat
+                frobeniusPower(i) := reduce lift frobeniusPower(i)
+          frobeniusPower := computeFrobeniusPowers()
+        m
+      modulus == m
+      if R has Finite then
+         size == d * size$R
+         random == UnVectorise([random()$R for i in 0..d1])
+      0 == 0$Rep
+      1 == 1$Rep
+      c * x == c *$Rep x
+      n * x == (n::R) *$Rep x
+      coerce(c:R):% == monomial(c,0)$Rep
+      coerce(x:%):OutputForm == coerce(x)$Rep
+      coefficient(x,e):R == coefficient(x,e)$Rep
+      reductum(x) == reductum(x)$Rep
+      leadingCoefficient x == (leadingCoefficient x)$Rep
+      degree x == (degree x)$Rep
+      lift(x) == x pretend Rep
+      reduce(p) == (monicDivide(p,m)$Rep).remainder
+      coerce(p) == reduce(p)
+      x = y == x =$Rep y
+      x + y == x +$Rep y
+      - x == -$Rep x
+      x * y ==
+        p := x *$Rep y
+        ans:=0$Rep
+        while (n:=degree p)>d1 repeat
+           ans:=ans + leadingCoefficient(p)*power.(n-d)
+           p := reductum p
+        ans+p
+      Vectorise(x) == [coefficient(lift(x),i) for i in 0..d1]
+      UnVectorise(vect) ==
+        reduce(+/[monomial(vect.(i+1),i) for i in 0..d1])
+      computePowers ==
+           mat : PrimitiveArray(%):= new(d,0)
+           mat.0:= reductum(-m)$Rep
+           w: % := monomial$Rep (1,1)
+           for i in 1..d1 repeat
+              mat.i := w *$Rep mat.(i-1)
+              if degree mat.i=d then
+                mat.i:= reductum mat.i + leadingCoefficient mat.i * mat.0
+           mat
+      if frobenius? then
+          computeFrobeniusPowers() ==
+            mat : PrimitiveArray(%):= new(d,1)
+            mat.1:= mult := monomial(1, size$R)$%
+            for i in 2..d1 repeat
+               mat.i := mult * mat.(i-1)
+            mat
+
+          frobenius(a:%):% ==
+            aq:% := 0
+            while a^=0 repeat
+              aq:= aq + leadingCoefficient(a)*frobeniusPower(degree a)
+              a := reductum a
+            aq
+         
+      pow == power
+      monomial(c,e)==
+         if e<d then monomial(c,e)$Rep
+         else
+            if e<=twod then
+               c * power.(e-d)
+            else
+               k1:=e quo twod
+               k2 := (e-k1*twod)::NonNegativeInteger
+               reduce((power.d1 **k1)*monomial(c,k2))
+      if R has Field then
+
+         (x:% exquo y:%):Union(%, "failed") ==
+            uv := extendedEuclidean(y, modulus(), x)$Rep
+            uv case "failed" => "failed"
+            return reduce(uv.coef1)
+
+         recip(y:%):Union(%, "failed") ==  1 exquo y
+         divide(x:%, y:%) ==
+            (q := (x exquo y)) case "failed" => error "not divisible"
+            [q, 0]
+
+--     An(MM) == Vectorise(-(reduce(reductum(m))::MM))
+--     LinearTransf(vect,MM) ==
+--       ans:= 0::SquareMatrix<d>(R)
+--       for i in 1..d do setelt(ans,i,1,vect.i)
+--       for j in 2..d do
+--          setelt(ans,1,j, elt(ans,d,j-1) * An(MM).1)
+--          for i in 2..d do
+--            setelt(ans,i,j, elt(ans,i-1,j-1) + elt(ans,d,j-1) * An(MM).i)
+--       ans
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MODFIELD ModularField}
+\pagehead{ModularField}{MODFIELD}
+\pagepic{ps/v103modularfield.ps}{MODFIELD}{1.00}
+See also:\\
+\refto{ModularRing}{MODRING}
+\refto{EuclideanModularRing}{EMR}
+<<domain MODFIELD ModularField>>=
+)abbrev domain MODFIELD ModularField
+++ These domains are used for the factorization and gcds
+++ of univariate polynomials over the integers in order to work modulo
+++ different  primes.
+++ See \spadtype{ModularRing}, \spadtype{EuclideanModularRing} 
+ModularField(R,Mod,reduction:(R,Mod) -> R,
+               merge:(Mod,Mod) -> Union(Mod,"failed"),
+                      exactQuo : (R,R,Mod) -> Union(R,"failed")) : C == T
+ where
+  R    :  CommutativeRing
+  Mod  :  AbelianMonoid
+
+  C == Field with
+                modulus :   %     -> Mod
+                       ++ modulus(x) \undocumented
+                coerce  :   %     -> R
+                       ++ coerce(x) \undocumented
+                reduce  : (R,Mod) -> %
+                       ++ reduce(r,m) \undocumented
+                exQuo   :  (%,%)  -> Union(%,"failed")
+                       ++ exQuo(x,y) \undocumented
+
+  T == ModularRing(R,Mod,reduction,merge,exactQuo)
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MODRING ModularRing}
+\pagehead{ModularRing}{MODRING}
+\pagepic{ps/v103modularring.ps}{MODRING}{1.00}
+See also:\\
+\refto{EuclideanModularRing}{EMR}
+\refto{ModularField}{MODFIELD}
+<<domain MODRING ModularRing>>=
+)abbrev domain MODRING ModularRing
+++ Author: P.Gianni, B.Trager
+++ Date Created:
+++ Date Last Updated:
+++ Basic Functions:
+++ Related Constructors:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ References:
+++ Description:
+++ These domains are used for the factorization and gcds
+++ of univariate polynomials over the integers in order to work modulo
+++ different  primes.
+++ See \spadtype{EuclideanModularRing} ,\spadtype{ModularField}
+
+ModularRing(R,Mod,reduction:(R,Mod) -> R,
+               merge:(Mod,Mod) -> Union(Mod,"failed"),
+                      exactQuo : (R,R,Mod) -> Union(R,"failed")) : C == T
+ where
+  R    :  CommutativeRing
+  Mod  :  AbelianMonoid
+
+  C == Ring with
+                modulus :   %     -> Mod
+                       ++ modulus(x) \undocumented
+                coerce  :   %     -> R
+                       ++ coerce(x) \undocumented
+                reduce  : (R,Mod) -> %
+                       ++ reduce(r,m) \undocumented
+                exQuo   :  (%,%)  -> Union(%,"failed")
+                       ++ exQuo(x,y) \undocumented
+                recip   :    %    -> Union(%,"failed")
+                       ++ recip(x) \undocumented
+                inv     :    %    -> %
+                       ++ inv(x) \undocumented
+
+  T == add
+    --representation
+      Rep:= Record(val:R,modulo:Mod)
+    --declarations
+      x,y: %
+
+    --define
+      modulus(x)   == x.modulo
+      coerce(x)    == x.val
+      coerce(i:Integer):% == [i::R,0]$Rep
+      i:Integer * x:% == (i::%)*x
+      coerce(x):OutputForm == (x.val)::OutputForm
+      reduce (a:R,m:Mod) == [reduction(a,m),m]$Rep
+
+      characteristic():NonNegativeInteger == characteristic()$R
+      0 == [0$R,0$Mod]$Rep
+      1 == [1$R,0$Mod]$Rep
+      zero? x == zero? x.val
+--      one? x == one? x.val
+      one? x == (x.val = 1)
+
+      newmodulo(m1:Mod,m2:Mod) : Mod ==
+        r:=merge(m1,m2)
+        r case "failed" => error "incompatible moduli"
+        r::Mod
+
+      x=y ==
+        x.val = y.val => true
+        x.modulo = y.modulo => false
+        (x-y).val = 0
+      x+y == reduce((x.val +$R y.val),newmodulo(x.modulo,y.modulo))
+      x-y == reduce((x.val -$R y.val),newmodulo(x.modulo,y.modulo))
+      -x  == reduce ((-$R x.val),x.modulo)
+      x*y == reduce((x.val *$R y.val),newmodulo(x.modulo,y.modulo))
+
+      exQuo(x,y) ==
+        xm:=x.modulo
+        if xm ^=$Mod y.modulo then xm:=newmodulo(xm,y.modulo)
+        r:=exactQuo(x.val,y.val,xm)
+        r case "failed"=> "failed"
+        [r::R,xm]$Rep
+
+      --if R has EuclideanDomain then
+      recip x ==
+        r:=exactQuo(1$R,x.val,x.modulo)
+        r case "failed" => "failed"
+        [r,x.modulo]$Rep
+
+      inv x ==
+        if (u:=recip x) case "failed" then error("not invertible")
+        else u::%
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MODMONOM ModuleMonomial}
+\pagehead{ModuleMonomial}{MODMONOM}
+\pagepic{ps/v103modulemonomial.ps}{MODMONOM}{1.00}
+See also:\\
+\refto{GeneralModulePolynomial}{GMODPOL}
+<<domain MODMONOM ModuleMonomial>>=
+)abbrev domain MODMONOM ModuleMonomial
+++ Description:
+++ This package \undocumented
+ModuleMonomial(IS: OrderedSet,
+               E: SetCategory,
+               ff:(MM, MM) -> Boolean): T == C where
+
+   MM ==> Record(index:IS, exponent:E)
+
+   T == OrderedSet  with
+        exponent: $ -> E
+               ++ exponent(x) \undocumented
+        index: $ -> IS
+               ++ index(x) \undocumented
+        coerce: MM -> $
+               ++ coerce(x) \undocumented
+        coerce: $ -> MM
+               ++ coerce(x) \undocumented
+        construct: (IS, E) -> $
+               ++ construct(i,e) \undocumented
+   C == MM  add
+        Rep:= MM
+        x:$ < y:$ == ff(x::Rep, y::Rep)
+        exponent(x:$):E == x.exponent
+        index(x:$): IS == x.index
+        coerce(x:$):MM == x::Rep::MM
+        coerce(x:MM):$ == x::Rep::$
+        construct(i:IS, e:E):$ == [i, e]$MM::Rep::$
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MOEBIUS MoebiusTransform}
+\pagehead{MoebiusTransform}{MOEBIUS}
+\pagepic{ps/v103moebiustransform.ps}{MOEBIUS}{1.00}
+<<domain MOEBIUS MoebiusTransform>>=
+)abbrev domain MOEBIUS MoebiusTransform
+++ 2-by-2 matrices acting on P1(F).
+++ Author: Stephen "Say" Watt
+++ Date Created: January 1987
+++ Date Last Updated: 11 April 1990
+++ Keywords:
+++ Examples:
+++ References:
+MoebiusTransform(F): Exports == Implementation where
+  ++ MoebiusTransform(F) is the domain of fractional linear (Moebius)
+  ++ transformations over F.
+  F : Field
+  OUT ==> OutputForm
+  P1F ==> OnePointCompletion F         -- projective 1-space over F
+ 
+  Exports ==> Group with
+ 
+    moebius: (F,F,F,F) -> %
+      ++ moebius(a,b,c,d) returns \spad{matrix [[a,b],[c,d]]}.
+    shift: F -> %
+      ++ shift(k) returns \spad{matrix [[1,k],[0,1]]} representing the map 
+      ++ \spad{x -> x + k}.
+    scale: F -> %
+      ++ scale(k) returns \spad{matrix [[k,0],[0,1]]} representing the map 
+      ++ \spad{x -> k * x}.
+    recip: () -> %
+      ++ recip() returns \spad{matrix [[0,1],[1,0]]} representing the map 
+      ++ \spad{x -> 1 / x}.
+    shift: (%,F) -> %
+      ++ shift(m,h) returns \spad{shift(h) * m} 
+      ++ (see \spadfunFrom{shift}{MoebiusTransform}).
+    scale: (%,F) -> %
+      ++ scale(m,h) returns \spad{scale(h) * m}
+      ++ (see \spadfunFrom{shift}{MoebiusTransform}).
+    recip: % -> %
+      ++ recip(m) = recip() * m
+    eval: (%,F) -> F
+      ++ eval(m,x) returns \spad{(a*x + b)/(c*x + d)} 
+      ++ where \spad{m = moebius(a,b,c,d)}
+      ++ (see \spadfunFrom{moebius}{MoebiusTransform}).
+    eval: (%,P1F) -> P1F
+      ++ eval(m,x) returns \spad{(a*x + b)/(c*x + d)} 
+      ++ where \spad{m = moebius(a,b,c,d)}
+      ++ (see \spadfunFrom{moebius}{MoebiusTransform}).
+
+  Implementation ==> add
+ 
+    Rep := Record(a: F,b: F,c: F,d: F)
+ 
+    moebius(aa,bb,cc,dd) == [aa,bb,cc,dd]
+ 
+    a(t:%):F == t.a
+    b(t:%):F == t.b
+    c(t:%):F == t.c
+    d(t:%):F == t.d
+ 
+    1 == moebius(1,0,0,1)
+    t * s ==
+      moebius(b(t)*c(s) + a(t)*a(s), b(t)*d(s) + a(t)*b(s), _
+              d(t)*c(s) + c(t)*a(s), d(t)*d(s) + c(t)*b(s))
+    inv t == moebius(d(t),-b(t),-c(t),a(t))
+ 
+    shift f == moebius(1,f,0,1)
+    scale f == moebius(f,0,0,1)
+    recip() == moebius(0,1,1,0)
+ 
+    shift(t,f) == moebius(a(t) + f*c(t), b(t) + f*d(t), c(t), d(t))
+    scale(t,f) == moebius(f*a(t),f*b(t),c(t),d(t))
+    recip t    == moebius(c(t),d(t),a(t),b(t))
+ 
+    eval(t:%,f:F) == (a(t)*f + b(t))/(c(t)*f + d(t))
+    eval(t:%,f:P1F) ==
+      (ff := retractIfCan(f)@Union(F,"failed")) case "failed" =>
+        (a(t)/c(t)) :: P1F
+      zero?(den := c(t) * (fff := ff :: F) + d(t)) => infinity()
+      ((a(t) * fff + b(t))/den) :: P1F
+ 
+    coerce t ==
+      var := "%x" :: OUT
+      num := (a(t) :: OUT) * var + (b(t) :: OUT)
+      den := (c(t) :: OUT) * var + (d(t) :: OUT)
+      rarrow(var,num/den)
+ 
+    proportional?: (List F,List F) -> Boolean
+    proportional?(list1,list2) ==
+      empty? list1 => empty? list2
+      empty? list2 => false
+      zero? (x1 := first list1) =>
+        (zero? first list2) and proportional?(rest list1,rest list2)
+      zero? (x2 := first list2) => false
+      map(#1 / x1,list1) = map(#1 / x2,list2)
+ 
+    t = s ==
+      list1 : List F := [a(t),b(t),c(t),d(t)]
+      list2 : List F := [a(s),b(s),c(s),d(s)]
+      proportional?(list1,list2)
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MRING MonoidRing}
+\pagehead{MonoidRing}{MRING}
+\pagepic{ps/v103monoidring.ps}{MRING}{1.00}
+<<domain MRING MonoidRing>>=
+)abbrev domain MRING MonoidRing
+++ Authors: Stephan M. Watt; revised by Johannes Grabmeier
+++ Date Created: January 1986
+++ Date Last Updated: 14 December 1995, Mike Dewar
+++ Basic Operations: *, +, monomials, coefficients
+++ Related Constructors: Polynomial
+++ Also See:
+++ AMS Classifications:
+++ Keywords: monoid ring, group ring, polynomials in non-commuting
+++  indeterminates
+++ References:
+++ Description:
+++  \spadtype{MonoidRing}(R,M), implements the algebra
+++  of all maps from the monoid M to the commutative ring R with
+++  finite support.
+++  Multiplication of two maps f and g is defined
+++  to map an element c of M to the (convolution) sum over {\em f(a)g(b)}
+++  such that {\em ab = c}. Thus M can be identified with a canonical
+++  basis and the maps can also be considered as formal linear combinations
+++  of the elements in M. Scalar multiples of a basis element are called
+++  monomials. A prominent example is the class of polynomials
+++  where the monoid is a direct product of the natural numbers
+++  with pointwise addition. When M is
+++  \spadtype{FreeMonoid Symbol}, one gets polynomials
+++  in infinitely many non-commuting variables. Another application
+++  area is representation theory of finite groups G, where modules
+++  over \spadtype{MonoidRing}(R,G) are studied.
+
+MonoidRing(R: Ring, M: Monoid): MRcategory == MRdefinition where
+    Term ==> Record(coef: R, monom: M)
+
+    MRcategory ==> Join(Ring, RetractableTo M, RetractableTo R) with
+        monomial         : (R, M) -> %
+          ++ monomial(r,m) creates a scalar multiple of the basis element m.
+        coefficient : (%, M) -> R
+          ++ coefficient(f,m) extracts the coefficient of m in f with respect
+          ++ to the canonical basis M.
+        coerce:   List Term -> %
+          ++ coerce(lt) converts a list of terms and coefficients to a member 
of the domain.
+        terms       : % -> List Term
+          ++ terms(f) gives the list of non-zero coefficients combined
+          ++ with their corresponding basis element as records.
+          ++ This is the internal representation.
+        map         : (R -> R, %) -> %
+          ++ map(fn,u) maps function fn onto the coefficients
+          ++ of the non-zero monomials of u.
+        monomial?   : % -> Boolean
+          ++ monomial?(f) tests if f is a single monomial.
+        coefficients: % -> List R
+          ++ coefficients(f) lists all non-zero coefficients.
+        monomials: % -> List %
+           ++ monomials(f) gives the list of all monomials whose
+           ++ sum is f.
+        numberOfMonomials: % -> NonNegativeInteger
+           ++ numberOfMonomials(f) is the number of non-zero coefficients
+           ++ with respect to the canonical basis.
+        if R has CharacteristicZero then CharacteristicZero
+        if R has CharacteristicNonZero then CharacteristicNonZero
+        if R has CommutativeRing then Algebra(R)
+        if (R has Finite and M has Finite) then Finite
+        if M has OrderedSet then
+          leadingMonomial   : % -> M
+            ++ leadingMonomial(f) gives the monomial of f whose
+            ++ corresponding monoid element is the greatest
+            ++ among all those with non-zero coefficients.
+          leadingCoefficient: % -> R
+            ++ leadingCoefficient(f) gives the coefficient of f, whose
+            ++ corresponding monoid element is the greatest
+            ++ among all those with non-zero coefficients.
+          reductum          : % -> %
+            ++ reductum(f) is f minus its leading monomial.
+
+    MRdefinition ==> add
+        Ex ==> OutputForm
+        Cf ==> coef
+        Mn ==> monom
+
+        Rep  := List Term
+
+        coerce(x: List Term): % == x :: %
+
+        monomial(r:R, m:M)  ==
+          r = 0 => empty()
+          [[r, m]]
+
+        if (R has Finite and M has Finite) then
+          size() == size()$R ** size()$M
+
+          index k ==
+            -- use p-adic decomposition of k
+            -- coefficient of p**j determines coefficient of index(i+p)$M
+            i:Integer := k rem size()
+            p:Integer := size()$R
+            n:Integer := size()$M
+            ans:% := 0
+            for j in 0.. while i > 0 repeat
+              h := i rem p
+              -- we use index(p) = 0$R
+              if h ^= 0 then
+                c : R := index(h :: PositiveInteger)$R
+                m : M := index((j+n) :: PositiveInteger)$M
+                --ans := ans + c *$% m
+                ans := ans + monomial(c, m)$%
+              i := i quo p
+            ans
+
+          lookup(z : %) : PositiveInteger ==
+            -- could be improved, if M has OrderedSet
+            -- z = index lookup z, n = lookup index n
+            -- use p-adic decomposition of k
+            -- coefficient of p**j determines coefficient of index(i+p)$M
+            zero?(z) => size()$% pretend PositiveInteger
+            liTe : List Term := terms z  -- all non-zero coefficients
+            p  : Integer := size()$R
+            n  : Integer := size()$M
+            res : Integer := 0
+            for te in liTe repeat
+              -- assume that lookup(p)$R = 0
+              l:NonNegativeInteger:=lookup(te.Mn)$M
+              ex : NonNegativeInteger := (n=l => 0;l)
+              co : Integer := lookup(te.Cf)$R
+              res := res + co * p ** ex
+            res pretend PositiveInteger
+
+          random() == index( (1+(random()$Integer rem size()$%) )_
+            pretend PositiveInteger)$%
+
+        0                   == empty()
+        1                   == [[1, 1]]
+        terms a             == (copy a) pretend List(Term)
+        monomials a         == [[t] for t in a]
+        coefficients a      == [t.Cf for t in a]
+        coerce(m:M):%       == [[1, m]]
+        coerce(r:R): % ==
+        -- coerce of ring
+          r = 0 => 0
+          [[r,    1]]
+        coerce(n:Integer): % ==
+        -- coerce of integers
+          n = 0 => 0
+          [[n::R, 1]]
+        - a                 == [[ -t.Cf, t.Mn] for t in a]
+        if R has noZeroDivisors
+           then
+            (r:R) * (a:%) ==
+              r = 0 => 0
+              [[r*t.Cf, t.Mn] for t in a]
+           else
+            (r:R) * (a:%) ==
+              r = 0 => 0
+              [[rt, t.Mn] for t in a | (rt:=r*t.Cf) ^= 0]
+        if R has noZeroDivisors
+           then
+            (n:Integer) * (a:%) ==
+              n = 0 => 0
+              [[n*t.Cf, t.Mn] for t in a]
+           else
+            (n:Integer) * (a:%) ==
+              n = 0 => 0
+              [[nt, t.Mn] for t in a | (nt:=n*t.Cf) ^= 0]
+        map(f, a)           == [[ft, t.Mn] for t in a | (ft:=f(t.Cf)) ^= 0]
+        numberOfMonomials a == #a
+
+        retractIfCan(a:%):Union(M, "failed") ==
+--          one?(#a) and one?(a.first.Cf) => a.first.Mn
+          ((#a) = 1) and ((a.first.Cf) = 1) => a.first.Mn
+          "failed"
+
+        retractIfCan(a:%):Union(R, "failed") ==
+--          one?(#a) and one?(a.first.Mn) => a.first.Cf
+          ((#a) = 1) and ((a.first.Mn) = 1) => a.first.Cf
+          "failed"
+
+        if R has noZeroDivisors then
+          if M has Group then
+            recip a ==
+              lt := terms a
+              #lt ^= 1 => "failed"
+              (u := recip lt.first.Cf) case "failed" => "failed"
+              --(u::R) * inv lt.first.Mn
+              monomial((u::R), inv lt.first.Mn)$%
+          else
+            recip a ==
+              #a ^= 1 or a.first.Mn ^= 1 => "failed"
+              (u := recip a.first.Cf) case "failed" => "failed"
+              u::R::%
+
+        mkTerm(r:R, m:M):Ex ==
+            r=1 => m::Ex
+            r=0 or m=1 => r::Ex
+            r::Ex * m::Ex
+
+        coerce(a:%):Ex ==
+            empty? a => (0$Integer)::Ex
+            empty? rest a => mkTerm(a.first.Cf, a.first.Mn)
+            reduce(_+, [mkTerm(t.Cf, t.Mn) for t in a])$List(Ex)
+
+        if M has OrderedSet then -- we mean totally ordered
+            -- Terms are stored in decending order.
+            leadingCoefficient a == (empty? a => 0; a.first.Cf)
+            leadingMonomial a    == (empty? a => 1; a.first.Mn)
+            reductum a           == (empty? a => a; rest a)
+
+            a = b ==
+                #a ^= #b => false
+                for ta in a for tb in b repeat
+                    ta.Cf ^= tb.Cf or ta.Mn ^= tb.Mn => return false
+                true
+
+            a + b ==
+                c:% := empty()
+                while not empty? a and not empty? b repeat
+                  ta := first a; tb := first b
+                  ra := rest a;  rb := rest b
+                  c :=
+                    ta.Mn > tb.Mn => (a := ra; concat_!(c, ta))
+                    ta.Mn < tb.Mn => (b := rb; concat_!(c, tb))
+                    a := ra; b := rb
+                    not zero?(r := ta.Cf+tb.Cf) =>
+                                        concat_!(c, [r, ta.Mn])
+                    c
+                concat_!(c, concat(a, b))
+
+            coefficient(a, m) ==
+                for t in a repeat
+                    if t.Mn = m then return t.Cf
+                    if t.Mn < m then return 0
+                0
+
+
+            if M has OrderedMonoid then
+
+            -- we use that multiplying an ordered list of monoid elements
+            -- by a single element respects the ordering
+
+              if R has noZeroDivisors then
+                a:% * b:% ==
+                  +/[[[ta.Cf*tb.Cf, ta.Mn*tb.Mn]$Term
+                    for tb in b ] for ta in reverse a]
+              else
+                a:% * b:% ==
+                  +/[[[r, ta.Mn*tb.Mn]$Term
+                    for tb in b | not zero?(r := ta.Cf*tb.Cf)]
+                      for ta in reverse a]
+            else -- M hasn't OrderedMonoid
+
+            -- we cannot assume that mutiplying an ordered list of
+            -- monoid elements by a single element respects the ordering:
+            -- we have to order and to collect equal terms
+              ge : (Term,Term) -> Boolean
+              ge(s,t) == t.Mn <= s.Mn
+
+              sortAndAdd : List Term -> List Term
+              sortAndAdd(liTe) ==  -- assume liTe not empty
+                liTe := sort(ge,liTe)
+                m : M :=  (first liTe).Mn
+                cf : R := (first liTe).Cf
+                res : List Term := []
+                for te in rest liTe repeat
+                  if m = te.Mn then
+                    cf := cf + te.Cf
+                  else
+                    if not zero? cf then res := cons([cf,m]$Term, res)
+                    m := te.Mn
+                    cf := te.Cf
+                if not zero? cf then res := cons([cf,m]$Term, res)
+                reverse res
+
+
+              if R has noZeroDivisors then
+                a:% * b:% ==
+                  zero? a => a
+                  zero? b => b  -- avoid calling sortAndAdd with []
+                  +/[sortAndAdd [[ta.Cf*tb.Cf, ta.Mn*tb.Mn]$Term
+                    for tb in b ] for ta in reverse a]
+              else
+                a:% * b:% ==
+                  zero? a => a
+                  zero? b => b  -- avoid calling sortAndAdd with []
+                  +/[sortAndAdd [[r, ta.Mn*tb.Mn]$Term
+                    for tb in b | not zero?(r := ta.Cf*tb.Cf)]
+                      for ta in reverse a]
+
+
+        else -- M hasn't OrderedSet
+            -- Terms are stored in random order.
+          a = b ==
+            #a ^= #b => false
+            brace(a pretend List(Term)) =$Set(Term) brace(b pretend List(Term))
+
+          coefficient(a, m) ==
+            for t in a repeat
+              t.Mn = m => return t.Cf
+            0
+
+          addterm(Tabl: AssociationList(M,R), r:R, m:M):R ==
+              (u := search(m, Tabl)) case "failed" => Tabl.m := r
+              zero?(r := r + u::R) => (remove_!(m, Tabl); 0)
+              Tabl.m := r
+
+          a + b ==
+              Tabl := table()$AssociationList(M,R)
+              for t in a repeat
+                  Tabl t.Mn := t.Cf
+              for t in b repeat
+                  addterm(Tabl, t.Cf, t.Mn)
+              [[Tabl m, m]$Term for m in keys Tabl]
+
+          a:% * b:% ==
+              Tabl := table()$AssociationList(M,R)
+              for ta in a repeat
+                  for tb in (b pretend List(Term)) repeat
+                      addterm(Tabl, ta.Cf*tb.Cf, ta.Mn*tb.Mn)
+              [[Tabl.m, m]$Term for m in keys Tabl]
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MSET Multiset}
+<<Multiset.input>>=
+-- mset.spad.pamphlet Multiset.input
+)spool Multiset.output
+)set message test on
+)set message auto off
+)clear all
+--S 1 of 14
+s := multiset [1,2,3,4,5,4,3,2,3,4,5,6,7,4,10]
+--R 
+--R
+--R   (1)  {1,2: 2,3: 3,4: 4,2: 5,6,7,10}
+--R                                               Type: Multiset 
PositiveInteger
+--E 1
+
+--S 2 of 14
+insert!(3,s)
+--R 
+--R
+--R   (2)  {1,2: 2,4: 3,4: 4,2: 5,6,7,10}
+--R                                               Type: Multiset 
PositiveInteger
+--E 2
+
+--S 3 of 14
+remove!(3,s,1)
+--R 
+--R
+--R   (3)  {1,2: 2,3: 3,4: 4,2: 5,6,7,10}
+--R                                               Type: Multiset 
PositiveInteger
+--E 3
+
+--S 4 of 14
+s
+--R 
+--R
+--R   (4)  {1,2: 2,3: 3,4: 4,2: 5,6,7,10}
+--R                                               Type: Multiset 
PositiveInteger
+--E 4
+
+--S 5 of 14
+remove!(5,s)
+--R 
+--R
+--R   (5)  {1,2: 2,3: 3,4: 4,6,7,10}
+--R                                               Type: Multiset 
PositiveInteger
+--E 5
+
+--S 6 of 14
+s
+--R 
+--R
+--R   (6)  {1,2: 2,3: 3,4: 4,6,7,10}
+--R                                               Type: Multiset 
PositiveInteger
+--E 6
+
+--S 7 of 14
+count(5,s)
+--R 
+--R
+--R   (7)  0
+--R                                                     Type: 
NonNegativeInteger
+--E 7
+
+--S 8 of 14
+t := multiset [2,2,2,-9]
+--R 
+--R
+--R   (8)  {3: 2,- 9}
+--R                                                       Type: Multiset 
Integer
+--E 8
+
+--S 9 of 14
+U := union(s,t)
+--R 
+--R
+--R   (9)  {1,5: 2,3: 3,4: 4,6,7,10,- 9}
+--R                                                       Type: Multiset 
Integer
+--E 9
+
+--S 10 of 14
+I := intersect(s,t)
+--R 
+--R
+--R   (10)  {5: 2}
+--R                                                       Type: Multiset 
Integer
+--E 10
+
+--S 11 of 14
+difference(s,t)
+--R 
+--R
+--R   (11)  {1,3: 3,4: 4,6,7,10}
+--R                                                       Type: Multiset 
Integer
+--E 11
+
+--S 12 of 14
+S := symmetricDifference(s,t)
+--R 
+--R
+--R   (12)  {1,3: 3,4: 4,6,7,10,- 9}
+--R                                                       Type: Multiset 
Integer
+--E 12
+
+--S 13 of 14
+(U = union(S,I))@Boolean
+--R 
+--R
+--R   (13)  true
+--R                                                                Type: 
Boolean
+--E 13
+
+--S 14 of 14
+t1 := multiset [1,2,2,3]; [t1 < t, t1 < s, t < s, t1 <= s]
+--R 
+--R
+--R   (14)  [false,true,false,true]
+--R                                                           Type: List 
Boolean
+--E 14
+)spool
+)lisp (bye)
+@
+<<Multiset.help>>=
+====================================================================
+Multiset examples
+====================================================================
+
+The domain Multiset(R) is similar to Set(R) except that multiplicities
+(counts of duplications) are maintained and displayed.  Use the
+operation multiset to create multisets from lists.  All the standard
+operations from sets are available for multisets.  An element with
+multiplicity greater than one has the multiplicity displayed first,
+then a colon, and then the element.
+
+Create a multiset of integers.
+
+  s := multiset [1,2,3,4,5,4,3,2,3,4,5,6,7,4,10]
+    {1,2: 2,3: 3,4: 4,2: 5,6,7,10}
+                          Type: Multiset PositiveInteger
+
+The operation insert! adds an element to a multiset.
+
+  insert!(3,s)
+    {1,2: 2,4: 3,4: 4,2: 5,6,7,10}
+                          Type: Multiset PositiveInteger
+
+Use remove! to remove an element.  If a third argument is present, it
+specifies how many instances to remove. Otherwise all instances of the
+element are removed.  Display the resulting multiset.
+
+  remove!(3,s,1); s
+
+    {1,2: 2,3: 3,4: 4,2: 5,6,7,10}
+                          Type: Multiset PositiveInteger
+
+  remove!(5,s); s
+    {1,2: 2,3: 3,4: 4,6,7,10}
+                          Type: Multiset PositiveInteger
+
+The operation count returns the number of copies of a given value.
+
+  count(5,s)
+    0
+                          Type: NonNegativeInteger
+
+A second multiset.
+
+  t := multiset [2,2,2,-9]
+    {3: 2,- 9}
+                          Type: Multiset Integer
+
+The union of two multisets is additive.
+
+  U := union(s,t)
+    {1,5: 2,3: 3,4: 4,6,7,10,- 9}
+                          Type: Multiset Integer
+
+The intersect operation gives the elements that are in common, with
+additive multiplicity.
+
+  I := intersect(s,t)
+    {5: 2}
+                          Type: Multiset Integer
+
+The difference of s and t consists of the elements that s has but t
+does not.  Elements are regarded as indistinguishable, so that if s
+and t have any element in common, the difference does not contain that
+element.
+
+  difference(s,t)
+    {1,3: 3,4: 4,6,7,10}
+                          Type: Multiset Integer
+
+The symmetricDifference is the union of difference(s, t) and difference(t, s).
+
+  S := symmetricDifference(s,t)
+    {1,3: 3,4: 4,6,7,10,- 9}
+                          Type: Multiset Integer
+
+Check that the union of the symmetricDifference and the intersect
+equals the union of the elements.
+
+  (U = union(S,I))@Boolean
+    true
+                          Type: Boolean
+
+Check some inclusion relations.
+
+  t1 := multiset [1,2,2,3]; [t1 < t, t1 < s, t < s, t1 <= s]
+    [false,true,false,true]
+                          Type: List Boolean
+
+See Also:
+o )show Multiset
+o $AXIOM/doc/src/algebra/mset.spad.dvi
+
+@
+\pagehead{Multiset}{MSET}
+\pagepic{ps/v103multiset.ps}{MSET}{1.00}
+<<domain MSET Multiset>>=
+)abbrev domain MSET Multiset
+++ Author:Stephen M. Watt, William H. Burge, Richard D. Jenks, Frederic Lehobey
+++ Date Created:NK
+++ Date Last Updated: 14 June 1994
+++ Basic Operations:
+++ Related Domains:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ Examples:
+++ References:
+++ Description: A multiset is a set with multiplicities.
+Multiset(S: SetCategory): MultisetAggregate S with
+        finiteAggregate
+        shallowlyMutable
+        multiset: () -> %
+          ++ multiset()$D creates an empty multiset of domain D.
+        multiset: S -> %
+          ++ multiset(s) creates a multiset with singleton s.
+        multiset: List S -> %
+          ++ multiset(ls) creates a multiset with elements from \spad{ls}.
+        members: % -> List S
+          ++ members(ms) returns a list of the elements of \spad{ms}
+          ++ {\em without} their multiplicity. See also \spadfun{parts}.
+        remove: (S,%,Integer) -> %
+          ++ remove(x,ms,number) removes at most \spad{number} copies of
+          ++ element x if \spad{number} is positive, all of them if
+          ++ \spad{number} equals zero, and all but at most \spad{-number} if
+          ++ \spad{number} is negative.
+        remove: ( S -> Boolean ,%,Integer) -> %
+          ++ remove(p,ms,number) removes at most \spad{number} copies of
+          ++ elements x such that \spad{p(x)} is \spadfun{true}
+          ++ if \spad{number} is positive, all of them if
+          ++ \spad{number} equals zero, and all but at most \spad{-number} if
+          ++ \spad{number} is negative.
+        remove_!: (S,%,Integer) -> %
+          ++ remove!(x,ms,number) removes destructively at most \spad{number}
+          ++ copies of element x if \spad{number} is positive, all
+          ++ of them if \spad{number} equals zero, and all but at most
+          ++ \spad{-number} if \spad{number} is negative.
+        remove_!: ( S -> Boolean ,%,Integer) -> %
+          ++ remove!(p,ms,number) removes destructively at most \spad{number}
+          ++ copies of elements x such that \spad{p(x)} is
+          ++ \spadfun{true} if \spad{number} is positive, all of them if
+          ++ \spad{number} equals zero, and all but at most \spad{-number} if
+          ++ \spad{number} is negative.
+
+    == add
+
+        Tbl ==> Table(S, Integer)
+        tbl ==> table$Tbl
+        Rep := Record(count: Integer, table: Tbl)
+
+        n: Integer
+        ms, m1, m2: %
+        t,  t1, t2: Tbl
+        D ==> Record(entry: S, count: NonNegativeInteger)
+        K ==> Record(key: S, entry: Integer)
+
+        elt(t:Tbl, s:S):Integer ==
+          a := search(s,t)$Tbl
+          a case "failed" => 0
+          a::Integer
+
+        empty():% == [0,tbl()]
+        multiset():% == empty()
+        dictionary():% == empty()                      -- DictionaryOperations
+        set():% == empty()
+        brace():% == empty()
+
+        construct(l:List S):% ==
+            t := tbl()
+            n := 0
+            for e in l repeat
+              t.e := inc t.e
+              n := inc n
+            [n, t]
+        multiset(l:List S):% == construct l
+        bag(l:List S):% == construct l                 -- BagAggregate
+        dictionary(l:List S):% == construct l          -- DictionaryOperations
+        set(l:List S):% == construct l
+        brace(l:List S):% == construct l
+
+        multiset(s:S):% == construct [s]
+
+        if S has ConvertibleTo InputForm then
+          convert(ms:%):InputForm ==
+            convert [convert("multiset"::Symbol)@InputForm,
+             convert(parts ms)@InputForm]
+
+        members(ms:%):List S == keys ms.table
+
+        coerce(ms:%):OutputForm ==
+            l: List OutputForm := empty()
+            t := ms.table
+            colon := ": " :: OutputForm
+            for e in keys t repeat
+                ex := e::OutputForm
+                n := t.e
+                item :=
+                  n > 1 => hconcat [n :: OutputForm,colon, ex]
+                  ex
+                l := cons(item,l)
+            brace l
+
+        duplicates(ms:%):List D ==                     -- MultiDictionary
+          ld : List D := empty()
+          t := ms.table
+          for e in keys t | (n := t.e) > 1 repeat
+            ld := cons([e,n::NonNegativeInteger],ld)
+          ld
+
+        extract_!(ms:%):S ==                           -- BagAggregate
+          empty? ms => error "extract: Empty multiset"
+          ms.count := dec ms.count
+          t := ms.table
+          e := inspect(t).key
+          if (n := t.e) > 1 then t.e := dec n
+           else remove_!(e,t)
+          e
+
+        inspect(ms:%):S == inspect(ms.table).key       -- BagAggregate
+
+        insert_!(e:S,ms:%):% ==                                -- BagAggregate
+            ms.count   := inc ms.count
+            ms.table.e := inc ms.table.e
+            ms
+
+        member?(e:S,ms:%):Boolean == member?(e,keys ms.table)
+
+        empty?(ms:%):Boolean == ms.count = 0
+
+        #(ms:%):NonNegativeInteger == ms.count::NonNegativeInteger
+
+        count(e:S, ms:%):NonNegativeInteger == ms.table.e::NonNegativeInteger
+
+        remove_!(e:S, ms:%, max:Integer):% ==
+          zero? max => remove_!(e,ms)
+          t := ms.table
+          if member?(e, keys t) then
+            ((n := t.e) <= max) =>
+              remove_!(e,t)
+              ms.count := ms.count-n
+            max > 0 =>
+              t.e := n-max
+              ms.count := ms.count-max
+            (n := n+max) > 0 =>
+              t.e := -max
+              ms.count := ms.count-n
+          ms
+
+        remove_!(p: S -> Boolean, ms:%, max:Integer):% ==
+          zero? max => remove_!(p,ms)
+          t := ms.table
+          for e in keys t | p(e) repeat
+            ((n := t.e) <= max) =>
+              remove_!(e,t)
+              ms.count := ms.count-n
+            max > 0 =>
+              t.e := n-max
+              ms.count := ms.count-max
+            (n := n+max) > 0 =>
+              t.e := -max
+              ms.count := ms.count-n
+          ms
+
+        remove(e:S, ms:%, max:Integer):% == remove_!(e, copy ms, max)
+
+        remove(p: S -> Boolean,ms:%,max:Integer):% == remove_!(p, copy ms, max)
+
+        remove_!(e:S, ms:%):% ==                       -- DictionaryOperations
+          t := ms.table
+          if member?(e, keys t) then
+            ms.count := ms.count-t.e
+            remove_!(e, t)
+          ms
+
+        remove_!(p:S ->Boolean, ms:%):% ==             -- DictionaryOperations
+          t := ms.table
+          for e in keys t | p(e) repeat
+            ms.count := ms.count-t.e
+            remove_!(e, t)
+          ms
+
+       select_!(p: S -> Boolean, ms:%):% ==            -- DictionaryOperations
+          remove_!(not p(#1), ms)
+
+        removeDuplicates_!(ms:%):% ==                  -- MultiDictionary
+          t := ms.table
+          l := keys t
+          for e in l repeat t.e := 1
+          ms.count := #l
+          ms
+
+        insert_!(e:S,ms:%,more:NonNegativeInteger):% ==        -- 
MultiDictionary
+            ms.count   := ms.count+more
+            ms.table.e := ms.table.e+more
+            ms
+
+        map_!(f: S->S, ms:%):% ==                      -- HomogeneousAggregate
+          t := ms.table
+          t1 := tbl()
+          for e in keys t repeat
+            t1.f(e) := t.e
+            remove_!(e, t)
+          ms.table := t1
+          ms
+
+       map(f: S -> S, ms:%):% == map_!(f, copy ms)     -- HomogeneousAggregate
+
+        parts(m:%):List S ==
+          l := empty()$List(S)
+          t := m.table
+          for e in keys t repeat
+            for i in 1..t.e repeat
+              l := cons(e,l)
+          l
+
+        union(m1:%, m2:%):% ==
+            t := tbl()
+            t1:= m1.table
+            t2:= m2.table
+            for e in keys t1 repeat t.e := t1.e
+            for e in keys t2 repeat t.e := t2.e + t.e
+            [m1.count + m2.count, t]
+
+        intersect(m1:%, m2:%):% ==
+--          if #m1 > #m2 then intersect(m2, m1)
+            t := tbl()
+            t1:= m1.table
+            t2:= m2.table
+            n := 0
+            for e in keys t1 repeat
+              m := min(t1.e,t2.e)
+              m > 0 =>
+                m := t1.e + t2.e
+                t.e := m
+                n := n + m
+            [n, t]
+
+        difference(m1:%, m2:%):% ==
+            t := tbl()
+            t1:= m1.table
+            t2:= m2.table
+            n := 0
+            for e in keys t1 repeat
+              k1 := t1.e
+              k2 := t2.e
+              k1 > 0 and k2 = 0 =>
+                t.e := k1
+                n := n + k1
+            n = 0 => empty()
+            [n, t]
+
+        symmetricDifference(m1:%, m2:%):% ==
+            union(difference(m1,m2), difference(m2,m1))
+
+        m1 = m2 ==
+            m1.count ^= m2.count => false
+            t1 := m1.table
+            t2 := m2.table
+            for e in keys t1 repeat
+                t1.e ^= t2.e => return false
+            for e in keys t2 repeat
+                t1.e ^= t2.e => return false
+            true
+
+        m1 < m2 ==
+            m1.count >= m2.count => false
+            t1 := m1.table
+            t2 := m2.table
+            for e in keys t1 repeat
+                t1.e > t2.e => return false
+            m1.count < m2.count
+
+        subset?(m1:%, m2:%):Boolean ==
+            m1.count > m2.count => false
+            t1 := m1.table
+            t2 := m2.table
+            for e in keys t1 repeat t1.e > t2.e => return false
+            true
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain MPOLY MultivariatePolynomial}
+<<MultivariatePolynomial.input>>=
+-- multpoly.spad.pamphlet MultivariatePolynomial.input
+)spool MultivariatePolynomial.output
+)set message test on
+)set message auto off
+)clear all
+--S 1 of 9
+m : MPOLY([x,y],INT) := (x^2 - x*y^3 +3*y)^2
+--R 
+--R
+--R         4     3 3     6       2     4      2
+--R   (1)  x  - 2y x  + (y  + 6y)x  - 6y x + 9y
+--R                                  Type: 
MultivariatePolynomial([x,y],Integer)
+--E 1
+
+--S 2 of 9
+m :: MPOLY([y,x],INT)
+--R 
+--R
+--R         2 6       4     3 3     2     2     4
+--R   (2)  x y  - 6x y  - 2x y  + 9y  + 6x y + x
+--R                                  Type: 
MultivariatePolynomial([y,x],Integer)
+--E 2
+
+--S 3 of 9
+p : MPOLY([x,y],POLY INT)
+--R 
+--R                                                                   Type: 
Void
+--E 3
+
+--S 4 of 9
+p :: POLY INT
+--R 
+--R
+--R   (4)  p
+--R                                                     Type: Polynomial 
Integer
+--E 4
+
+--S 5 of 9
+% :: MPOLY([a,b],POLY INT)
+--R 
+--R
+--R   (5)  p
+--R                       Type: MultivariatePolynomial([a,b],Polynomial 
Integer)
+--E 5
+
+--S 6 of 9
+q : UP(x, FRAC MPOLY([y,z],INT))
+--R 
+--R                                                                   Type: 
Void
+--E 6
+
+--S 7 of 9
+q := (x^2 - x*(z+1)/y +2)^2 
+--R 
+--R
+--R                             2    2
+--R         4   - 2z - 2  3   4y  + z  + 2z + 1  2   - 4z - 4
+--R   (7)  x  + -------- x  + ----------------- x  + -------- x + 4
+--R                 y                  2                 y
+--R                                   y
+--R Type: UnivariatePolynomial(x,Fraction 
MultivariatePolynomial([y,z],Integer))
+--E 7
+
+--S 8 of 9
+q :: UP(z, FRAC MPOLY([x,y],INT))
+--R 
+--R
+--R   (8)
+--R    2            3     2             2 4       3      2      2            2
+--R   x   2   - 2y x  + 2x  - 4y x     y x  - 2y x  + (4y  + 1)x  - 4y x + 4y
+--R   -- z  + -------------------- z + ---------------------------------------
+--R    2                2                                  2
+--R   y                y                                  y
+--R Type: UnivariatePolynomial(z,Fraction 
MultivariatePolynomial([x,y],Integer))
+--E 8
+
+--S 9 of 9
+q :: MPOLY([x,z], FRAC UP(y,INT))
+--R 
+--R
+--R                                               2
+--R         4      2     2  3     1  2    2     4y  + 1  2      4     4
+--R   (9)  x  + (- - z - -)x  + (-- z  + -- z + -------)x  + (- - z - -)x + 4
+--R                y     y        2       2         2           y     y
+--R                              y       y         y
+--R Type: MultivariatePolynomial([x,z],Fraction 
UnivariatePolynomial(y,Integer))
+--E 9
+)spool
+)lisp (bye)
+@
+<<MultivariatePolynomial.help>>=
+====================================================================
+MultivariatePolynomial examples
+====================================================================
+
+The domain constructor MultivariatePolynomial is similar to Polynomial
+except that it specifies the variables to be used.  Polynomial are
+available for MultivariatePolynomial.  The abbreviation for
+MultivariatePolynomial is MPOLY.  The type expressions
+
+  MultivariatePolynomial([x,y],Integer)
+  MPOLY([x,y],INT) 
+
+refer to the domain of multivariate polynomials in the variables x and
+y where the coefficients are restricted to be integers.  The first
+variable specified is the main variable and the display of the polynomial
+reflects this.
+
+This polynomial appears with terms in descending powers of the variable x.
+
+  m : MPOLY([x,y],INT) := (x^2 - x*y^3 +3*y)^2
+     4     3 3     6       2     4      2
+    x  - 2y x  + (y  + 6y)x  - 6y x + 9y
+                     Type: MultivariatePolynomial([x,y],Integer)
+
+It is easy to see a different variable ordering by doing a conversion.
+
+  m :: MPOLY([y,x],INT)
+     2 6       4     3 3     2     2     4
+    x y  - 6x y  - 2x y  + 9y  + 6x y + x
+                     Type: MultivariatePolynomial([y,x],Integer)
+
+You can use other, unspecified variables, by using Polynomial in the
+coefficient type of MPOLY.
+
+  p : MPOLY([x,y],POLY INT)
+                     Type: Void
+
+Conversions can be used to re-express such polynomials in terms of
+the other variables.  For example, you can first push all the
+variables into a polynomial with integer coefficients.
+
+  p :: POLY INT
+    p
+                     Type: Polynomial Integer
+
+Now pull out the variables of interest.
+
+  % :: MPOLY([a,b],POLY INT)
+    p
+                     Type: MultivariatePolynomial([a,b],Polynomial Integer)
+
+Restriction:
+  Axiom does not allow you to create types where MultivariatePolynomial
+  is contained in the coefficient type of Polynomial. Therefore,
+  MPOLY([x,y],POLY INT) is legal but POLY MPOLY([x,y],INT) is not.
+
+Multivariate polynomials may be combined with univariate polynomials
+to create types with special structures.
+
+  q : UP(x, FRAC MPOLY([y,z],INT))
+                          Type: Void
+
+This is a polynomial in x whose coefficients are quotients of polynomials 
+in y and z.
+
+  q := (x^2 - x*(z+1)/y +2)^2 
+                             2    2
+         4   - 2z - 2  3   4y  + z  + 2z + 1  2   - 4z - 4
+   (7)  x  + -------- x  + ----------------- x  + -------- x + 4
+                 y                  2                 y
+                                   y
+ Type: UnivariatePolynomial(x,Fraction MultivariatePolynomial([y,z],Integer))
+
+Use conversions for structural rearrangements. z does not appear in a
+denominator and so it can be made the main variable.
+
+  q :: UP(z, FRAC MPOLY([x,y],INT))
+    2            3     2             2 4       3      2      2            2
+   x   2   - 2y x  + 2x  - 4y x     y x  - 2y x  + (4y  + 1)x  - 4y x + 4y
+   -- z  + -------------------- z + ---------------------------------------
+    2                2                                  2
+   y                y                                  y
+ Type: UnivariatePolynomial(z,Fraction MultivariatePolynomial([x,y],Integer))
+
+Or you can make a multivariate polynomial in x and z whose
+coefficients are fractions in polynomials in y.
+
+  q :: MPOLY([x,z], FRAC UP(y,INT))
+     4      2     2  3     1  2    2     4y  + 1  2      4     4
+    x  + (- - z - -)x  + (-- z  + -- z + -------)x  + (- - z - -)x + 4
+            y     y        2       2         2           y     y
+                          y       y         y
+ Type: MultivariatePolynomial([x,z],Fraction UnivariatePolynomial(y,Integer))
+
+A conversion like q :: MPOLY([x,y], FRAC UP(z,INT)) is not possible in
+this example because y appears in the denominator of a fraction.  As
+you can see, Axiom provides extraordinary flexibility in the
+manipulation and display of expressions via its conversion facility.
+
+See Also:
+o )help DistributedMultivariatePolynomial
+o )help UnivariatePolynomial
+o )help Polynomial
+o )show MultivariatePolynomial
+o $AXIOM/doc/src/algebra/multpoly.spad.dvi
+
+@
+\pagehead{MultivariatePolynomial}{MPOLY}
+\pagepic{ps/v103multivariatepolynomial.ps}{MPOLY}{1.00}
+See also:\\
+\refto{Polynomial}{POLY}
+\refto{SparseMultivariatePolynomial}{SMP}
+\refto{IndexedExponents}{INDE}
+<<domain MPOLY MultivariatePolynomial>>=
+)abbrev domain MPOLY MultivariatePolynomial
+++ Author: Dave Barton, Barry Trager
+++ Date Created:
+++ Date Last Updated:
+++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
+++ resultant, gcd
+++ Related Constructors: SparseMultivariatePolynomial, Polynomial
+++ Also See:
+++ AMS Classifications:
+++ Keywords: polynomial, multivariate
+++ References:
+++ Description:
+++   This type is the basic representation of sparse recursive multivariate
+++ polynomials whose variables are from a user specified list of symbols.
+++ The ordering is specified by the position of the variable in the list.
+++ The coefficient ring may be non commutative,
+++ but the variables are assumed to commute.
+
+MultivariatePolynomial(vl:List Symbol, R:Ring)
+   ==  SparseMultivariatePolynomial(--SparseUnivariatePolynomial,
+           R, OrderedVariableList vl)
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \chapter{Chapter N}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain NONE None}
@@ -38318,6 +41398,846 @@ PlaneAlgebraicCurvePlot(): 
PlottablePlaneCurveCategory _
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain POLY Polynomial}
+<<Polynomial.input>>=
+-- multpoly.spad.pamphlet Polynomial.input
+)spool Polynomial.output
+)set message test on
+)set message auto off
+--S 1 of 46
+x + 1
+--R 
+--R
+--R   (1)  x + 1
+--R                                                     Type: Polynomial 
Integer
+--E 1
+
+--S 2 of 46
+z - 2.3
+--R 
+--R
+--R   (2)  z - 2.3
+--R                                                       Type: Polynomial 
Float
+--E 2
+
+--S 3 of 46
+y**2 - z + 3/4
+--R 
+--R
+--R               2   3
+--R   (3)  - z + y  + -
+--R                   4
+--R                                            Type: Polynomial Fraction 
Integer
+--E 3
+
+--S 4 of 46
+y **2 + x*y + y
+--R 
+--R
+--R         2
+--R   (4)  y  + (x + 1)y
+--R                                                     Type: Polynomial 
Integer
+--E 4
+
+--S 5 of 46
+% :: DMP([y,x],INT)
+--R 
+--R
+--R         2
+--R   (5)  y  + y x + y
+--R                       Type: 
DistributedMultivariatePolynomial([y,x],Integer)
+--E 5
+
+--S 6 of 46
+p := (y-1)**2 * x * z
+--R 
+--R
+--R            2
+--R   (6)  (x y  - 2x y + x)z
+--R                                                     Type: Polynomial 
Integer
+--E 6
+
+--S 7 of 46
+q := (y-1) * x * (z+5)
+--R 
+--R
+--R   (7)  (x y - x)z + 5x y - 5x
+--R                                                     Type: Polynomial 
Integer
+--E 7
+
+--S 8 of 46
+factor(q)
+--R 
+--R
+--R   (8)  x(y - 1)(z + 5)
+--R                                            Type: Factored Polynomial 
Integer
+--E 8
+
+--S 9 of 46
+p - q**2
+--R 
+--R
+--R   (9)
+--R         2 2     2     2  2          2      2       2             2
+--R     (- x y  + 2x y - x )z  + ((- 10x  + x)y  + (20x  - 2x)y - 10x  + x)z
+--R   + 
+--R          2 2      2       2
+--R     - 25x y  + 50x y - 25x
+--R                                                     Type: Polynomial 
Integer
+--E 9
+
+--S 10 of 46
+gcd(p,q)
+--R 
+--R
+--R   (10)  x y - x
+--R                                                     Type: Polynomial 
Integer
+--E 10
+
+--S 11 of 46
+factor %
+--R 
+--R
+--R   (11)  x(y - 1)
+--R                                            Type: Factored Polynomial 
Integer
+--E 11
+
+--S 12 of 46
+lcm(p,q)
+--R 
+--R
+--R             2             2        2
+--R   (12)  (x y  - 2x y + x)z  + (5x y  - 10x y + 5x)z
+--R                                                     Type: Polynomial 
Integer
+--E 12
+
+--S 13 of 46
+content p
+--R 
+--R
+--R   (13)  1
+--R                                                        Type: 
PositiveInteger
+--E 13
+
+--S 14 of 46
+resultant(p,q,z)
+--R 
+--R
+--R           2 3      2 2      2      2
+--R   (14)  5x y  - 15x y  + 15x y - 5x
+--R                                                     Type: Polynomial 
Integer
+--E 14
+
+--S 15 of 46
+resultant(p,q,x)
+--R 
+--R
+--R   (15)  0
+--R                                                     Type: Polynomial 
Integer
+--E 15
+
+--S 16 of 46
+mainVariable p
+--R 
+--R
+--R   (16)  z
+--R                                                      Type: 
Union(Symbol,...)
+--E 16
+
+--S 17 of 46
+mainVariable(1 :: POLY INT)
+--R 
+--R
+--R   (17)  "failed"
+--R                                                    Type: 
Union("failed",...)
+--E 17
+
+--S 18 of 46
+ground? p
+--R 
+--R
+--R   (18)  false
+--R                                                                Type: 
Boolean
+--E 18
+
+--S 19 of 46
+ground?(1 :: POLY INT)
+--R 
+--R
+--R   (19)  true
+--R                                                                Type: 
Boolean
+--E 19
+
+--S 20 of 46
+variables p
+--R 
+--R
+--R   (20)  [z,y,x]
+--R                                                            Type: List 
Symbol
+--E 20
+
+--S 21 of 46
+degree(p,x)
+--R 
+--R
+--R   (21)  1
+--R                                                        Type: 
PositiveInteger
+--E 21
+
+--S 22 of 46
+degree(p,y)
+--R 
+--R
+--R   (22)  2
+--R                                                        Type: 
PositiveInteger
+--E 22
+
+--S 23 of 46
+degree(p,z)
+--R 
+--R
+--R   (23)  1
+--R                                                        Type: 
PositiveInteger
+--E 23
+
+--S 24 of 46
+degree(p,[x,y,z])
+--R 
+--R
+--R   (24)  [1,2,1]
+--R                                                Type: List 
NonNegativeInteger
+--E 24
+
+--S 25 of 46
+minimumDegree(p,z)
+--R 
+--R
+--R   (25)  1
+--R                                                        Type: 
PositiveInteger
+--E 25
+
+--S 26 of 46
+totalDegree p
+--R 
+--R
+--R   (26)  4
+--R                                                        Type: 
PositiveInteger
+--E 26
+
+--S 27 of 46
+leadingMonomial p
+--R 
+--R
+--R            2
+--R   (27)  x y z
+--R                                                     Type: Polynomial 
Integer
+--E 27
+
+--S 28 of 46
+reductum p
+--R 
+--R
+--R   (28)  (- 2x y + x)z
+--R                                                     Type: Polynomial 
Integer
+--E 28
+
+--S 29 of 46
+p - leadingMonomial p - reductum p
+--R 
+--R
+--R   (29)  0
+--R                                                     Type: Polynomial 
Integer
+--E 29
+
+--S 30 of 46
+leadingCoefficient p
+--R 
+--R
+--R   (30)  1
+--R                                                        Type: 
PositiveInteger
+--E 30
+
+--S 31 of 46
+p
+--R 
+--R
+--R             2
+--R   (31)  (x y  - 2x y + x)z
+--R                                                     Type: Polynomial 
Integer
+--E 31
+
+--S 32 of 46
+eval(p,x,w)
+--R 
+--R
+--R             2
+--R   (32)  (w y  - 2w y + w)z
+--R                                                     Type: Polynomial 
Integer
+--E 32
+
+--S 33 of 46
+eval(p,x,1)
+--R 
+--R
+--R           2
+--R   (33)  (y  - 2y + 1)z
+--R                                                     Type: Polynomial 
Integer
+--E 33
+
+--S 34 of 46
+eval(p,x,y^2 - 1)
+--R 
+--R
+--R           4     3
+--R   (34)  (y  - 2y  + 2y - 1)z
+--R                                                     Type: Polynomial 
Integer
+--E 34
+
+--S 35 of 46
+D(p,x)
+--R 
+--R
+--R           2
+--R   (35)  (y  - 2y + 1)z
+--R                                                     Type: Polynomial 
Integer
+--E 35
+
+--S 36 of 46
+D(p,y)
+--R 
+--R
+--R   (36)  (2x y - 2x)z
+--R                                                     Type: Polynomial 
Integer
+--E 36
+
+--S 37 of 46
+D(p,z)
+--R 
+--R
+--R            2
+--R   (37)  x y  - 2x y + x
+--R                                                     Type: Polynomial 
Integer
+--E 37
+
+--S 38 of 46
+integrate(p,y)
+--R 
+--R
+--R          1    3      2
+--R   (38)  (- x y  - x y  + x y)z
+--R          3
+--R                                            Type: Polynomial Fraction 
Integer
+--E 38
+
+--S 39 of 46
+qr := monicDivide(p,x+1,x)
+--R 
+--R
+--R                      2                           2
+--R   (39)  [quotient= (y  - 2y + 1)z,remainder= (- y  + 2y - 1)z]
+--R     Type: Record(quotient: Polynomial Integer,remainder: Polynomial 
Integer)
+--E 39
+
+--S 40 of 46
+qr.remainder
+--R 
+--R
+--R             2
+--R   (40)  (- y  + 2y - 1)z
+--R                                                     Type: Polynomial 
Integer
+--E 40
+
+--S 41 of 46
+p - ((x+1) * qr.quotient + qr.remainder)
+--R 
+--R
+--R   (41)  0
+--R                                                     Type: Polynomial 
Integer
+--E 41
+
+--S 42 of 46
+p/q
+--R 
+--R
+--R         (y - 1)z
+--R   (42)  --------
+--R           z + 5
+--R                                            Type: Fraction Polynomial 
Integer
+--E 42
+
+--S 43 of 46
+(2/3) * x**2 - y + 4/5 
+--R 
+--R
+--R               2  2   4
+--R   (43)  - y + - x  + -
+--R               3      5
+--R                                            Type: Polynomial Fraction 
Integer
+--E 43
+
+--S 44 of 46
+% :: FRAC POLY INT
+--R 
+--R
+--R                    2
+--R         - 15y + 10x  + 12
+--R   (44)  -----------------
+--R                 15
+--R                                            Type: Fraction Polynomial 
Integer
+--E 44
+
+--S 45 of 46
+% :: POLY FRAC INT
+--R 
+--R
+--R               2  2   4
+--R   (45)  - y + - x  + -
+--R               3      5
+--R                                            Type: Polynomial Fraction 
Integer
+--E 45
+
+--S 46 of 46
+map(numeric,%)
+--R 
+--R
+--R                                            2
+--R   (46)  - 1.0 y + 0.6666666666 6666666667 x  + 0.8
+--R                                                       Type: Polynomial 
Float
+--E 46
+)spool
+)lisp (bye)
+@
+<<Polynomial.help>>=
+====================================================================
+Polynomial examples
+====================================================================
+
+The domain constructor Polynomial (abbreviation: POLY) provides
+polynomials with an arbitrary number of unspecified variables.
+
+It is used to create the default polynomial domains in Axiom.  Here
+the coefficients are integers.
+
+  x + 1
+    x + 1
+                          Type: Polynomial Integer
+
+Here the coefficients have type Float.
+
+  z - 2.3
+    z - 2.3
+                           Type: Polynomial Float
+
+And here we have a polynomial in two variables with coefficients which
+have type Fraction Integer.
+
+  y**2 - z + 3/4
+           2   3
+    - z + y  + -
+               4
+                           Type: Polynomial Fraction Integer
+
+The representation of objects of domains created by Polynomial is that
+of recursive univariate polynomials. The term univariate means "one
+variable". The term multivariate means "possibly more than one variable".
+
+This recursive structure is sometimes obvious from the display of a polynomial.
+
+  y **2 + x*y + y
+     2
+    y  + (x + 1)y
+                           Type: Polynomial Integer
+
+In this example, you see that the polynomial is stored as a polynomial
+in y with coefficients that are polynomials in x with integer coefficients.  
+In fact, you really don't need to worry about the representation unless you 
+are working on an advanced application where it is critical. The polynomial 
+types created from DistributedMultivariatePolynomial and
+NewDistributedMultivariatePolynomial are stored and displayed in a
+non-recursive manner.
+
+You see a "flat" display of the above polynomial by converting to
+one of those types.
+
+  % :: DMP([y,x],INT)
+     2
+    y  + y x + y
+                   Type: DistributedMultivariatePolynomial([y,x],Integer)
+
+We will demonstrate many of the polynomial facilities by using two
+polynomials with integer coefficients.
+
+By default, the interpreter expands polynomial expressions, even if they
+are written in a factored format.
+
+  p := (y-1)**2 * x * z
+        2
+    (x y  - 2x y + x)z
+                   Type: Polynomial Integer
+
+See Factored to see how to create objects in factored form directly.
+
+  q := (y-1) * x * (z+5)
+    (x y - x)z + 5x y - 5x
+                   Type: Polynomial Integer
+
+The fully factored form can be recovered by using factor.
+
+  factor(q)
+    x(y - 1)(z + 5)
+                   Type: Factored Polynomial Integer
+
+This is the same name used for the operation to factor integers.  Such
+reuse of names is called overloading and makes it much easier to think
+of solving problems in general ways.  Axiom facilities for factoring
+polynomials created with Polynomial are currently restricted to the
+integer and rational number coefficient cases.
+
+The standard arithmetic operations are available for polynomials.
+
+  p - q**2
+         2 2     2     2  2          2      2       2             2
+     (- x y  + 2x y - x )z  + ((- 10x  + x)y  + (20x  - 2x)y - 10x  + x)z
+   + 
+          2 2      2       2
+     - 25x y  + 50x y - 25x
+                   Type: Polynomial Integer
+
+The operation gcd is used to compute the greatest common divisor of
+two polynomials.
+
+  gcd(p,q)
+    x y - x
+                   Type: Polynomial Integer
+
+In the case of p and q, the gcd is obvious from their definitions.  We
+factor the gcd to show this relationship better.
+
+  factor %
+    x(y - 1)
+                   Type: Factored Polynomial Integer
+
+The least common multiple is computed by using lcm.
+
+  lcm(p,q)
+        2             2        2
+    (x y  - 2x y + x)z  + (5x y  - 10x y + 5x)z
+                   Type: Polynomial Integer
+
+Use content to compute the greatest common divisor of the coefficients
+of the polynomial.
+
+  content p
+    1
+                   Type: PositiveInteger
+
+Many of the operations on polynomials require you to specify a variable.  
+For example, resultant requires you to give the variable in which the 
+polynomials should be expressed.
+
+This computes the resultant of the values of p and q, considering them 
+as polynomials in the variable z.  They do not share a root when thought 
+of as polynomials in z.
+
+  resultant(p,q,z)
+       2 3      2 2      2      2
+     5x y  - 15x y  + 15x y - 5x
+                   Type: Polynomial Integer
+
+This value is 0 because as polynomials in x the polynomials have a
+common root.
+
+  resultant(p,q,x)
+    0
+                   Type: Polynomial Integer
+
+The data type used for the variables created by Polynomial is Symbol.
+As mentioned above, the representation used by Polynomial is recursive
+and so there is a main variable for nonconstant polynomials.
+
+The operation mainVariable returns this variable.  The return type is 
+actually a union of Symbol and "failed".
+
+  mainVariable p
+    z
+                   Type: Union(Symbol,...)
+
+The latter branch of the union is be used if the polynomial has no
+variables, that is, is a constant.
+
+  mainVariable(1 :: POLY INT)
+    "failed"
+                   Type: Union("failed",...)
+
+You can also use the predicate ground? to test whether a polynomial is
+in fact a member of its ground ring.
+
+  ground? p
+    false
+                   Type: Boolean
+
+  ground?(1 :: POLY INT)
+    true
+                   Type: Boolean
+
+The complete list of variables actually used in a particular polynomial 
+is returned by variables.  For constant polynomials, this list is empty.
+
+  variables p
+    [z,y,x]
+                   Type: List Symbol
+
+The degree operation returns the degree of a polynomial in a specific variable.
+
+  degree(p,x)
+    1
+                   Type: PositiveInteger
+
+  degree(p,y)
+    2
+                   Type: PositiveInteger
+
+  degree(p,z)
+    1
+                   Type: PositiveInteger
+
+If you give a list of variables for the second argument, a list of the
+degrees in those variables is returned.
+
+  degree(p,[x,y,z])
+    [1,2,1]
+                   Type: List NonNegativeInteger
+
+The minimum degree of a variable in a polynomial is computed using 
+minimumDegree.
+
+  minimumDegree(p,z)
+    1
+                   Type: PositiveInteger
+
+The total degree of a polynomial is returned by totalDegree.
+
+  totalDegree p
+    4
+                   Type: PositiveInteger
+
+It is often convenient to think of a polynomial as a leading monomial plus
+the remaining terms.
+
+  leadingMonomial p
+        2
+     x y z
+                   Type: Polynomial Integer
+
+The reductum operation returns a polynomial consisting of the sum of
+the monomials after the first.
+
+  reductum p
+    (- 2x y + x)z
+                   Type: Polynomial Integer
+
+These have the obvious relationship that the original polynomial is
+equal to the leading monomial plus the reductum.
+
+  p - leadingMonomial p - reductum p
+    0
+                   Type: Polynomial Integer
+
+The value returned by leadingMonomial includes the coefficient of that term.  
+This is extracted by using leadingCoefficient on the original polynomial.
+
+  leadingCoefficient p
+    1
+                   Type: PositiveInteger
+
+The operation eval is used to substitute a value for a variable in a 
+polynomial.
+
+  p
+        2
+    (x y  - 2x y + x)z
+                   Type: Polynomial Integer
+
+This value may be another variable, a constant or a polynomial.
+
+  eval(p,x,w)
+        2
+    (w y  - 2w y + w)z
+                   Type: Polynomial Integer
+
+  eval(p,x,1)
+      2
+    (y  - 2y + 1)z
+                   Type: Polynomial Integer
+
+Actually, all the things being substituted are just polynomials,
+some more trivial than others.
+
+  eval(p,x,y^2 - 1)
+      4     3
+    (y  - 2y  + 2y - 1)z
+                   Type: Polynomial Integer
+
+Derivatives are computed using the D operation.
+
+  D(p,x)
+      2
+    (y  - 2y + 1)z
+                   Type: Polynomial Integer
+
+The first argument is the polynomial and the second is the variable.
+
+  D(p,y)
+    (2x y - 2x)z
+                   Type: Polynomial Integer
+
+Even if the polynomial has only one variable, you must specify it.
+
+  D(p,z)
+       2
+    x y  - 2x y + x
+                   Type: Polynomial Integer
+
+Integration of polynomials is similar and the integrate operation is used.
+
+Integration requires that the coefficients support division. Axiom
+converts polynomials over the integers to polynomials over the rational 
+numbers before integrating them.
+
+  integrate(p,y)
+     1    3      2
+    (- x y  - x y  + x y)z
+     3
+                  Type: Polynomial Fraction Integer
+
+It is not possible, in general, to divide two polynomials.  In our
+example using polynomials over the integers, the operation monicDivide
+divides a polynomial by a monic polynomial (that is, a polynomial with
+leading coefficient equal to 1).  The result is a record of the
+quotient and remainder of the division.
+
+You must specify the variable in which to express the polynomial.
+
+  qr := monicDivide(p,x+1,x)
+                 2                           2
+    [quotient= (y  - 2y + 1)z,remainder= (- y  + 2y - 1)z]
+     Type: Record(quotient: Polynomial Integer,remainder: Polynomial Integer)
+
+The selectors of the components of the record are quotient and remainder.  
+Issue this to extract the remainder.
+
+  qr.remainder
+        2
+    (- y  + 2y - 1)z
+                       Type: Polynomial Integer
+
+Now that we can extract the components, we can demonstrate the
+relationship among them and the arguments to our original expression
+qr := monicDivide(p,x+1,x).
+
+  p - ((x+1) * qr.quotient + qr.remainder)
+    0
+                       Type: Polynomial Integer
+
+If the / operator is used with polynomials, a fraction object is
+created.  In this example, the result is an object of type 
+Fraction Polynomial Integer.
+
+  p/q
+    (y - 1)z
+    --------
+      z + 5
+                       Type: Fraction Polynomial Integer
+
+If you use rational numbers as polynomial coefficients, the
+resulting object is of type Polynomial Fraction Integer.
+
+  (2/3) * x**2 - y + 4/5 
+          2  2   4
+    - y + - x  + -
+          3      5
+                       Type: Polynomial Fraction Integer
+
+This can be converted to a fraction of polynomials and back again, if
+required.
+
+  % :: FRAC POLY INT
+               2
+    - 15y + 10x  + 12
+    -----------------
+            15
+                       Type: Fraction Polynomial Integer
+
+  % :: POLY FRAC INT
+          2  2   4
+    - y + - x  + -
+          3      5
+                       Type: Polynomial Fraction Integer
+
+To convert the coefficients to floating point, map the numeric
+operation on the coefficients of the polynomial.
+
+  map(numeric,%)
+    - 1.0 y + 0.6666666666 6666666667 x  + 0.8
+                       Type: Polynomial Float
+
+See Also:
+o )help Factored
+o )help UnivariatePolynomial
+o )help MultivariatePolynomial
+o )help DistributedMultivariatePolynomial
+o )help NewDistributedMultivariatePolynomial
+o )show Polynomial
+o $AXIOM/doc/src/algebra/multpoly.spad.dvi
+
+@
+\pagehead{Polynomial}{POLY}
+\pagepic{ps/v103polynomial.ps}{POLY}{1.00}
+See also:\\
+\refto{MultivariatePolynomial}{MPOLY}
+\refto{SparseMultivariatePolynomial}{SMP}
+\refto{IndexedExponents}{INDE}
+<<domain POLY Polynomial>>=
+)abbrev domain POLY Polynomial
+++ Author: Dave Barton, Barry Trager
+++ Date Created:
+++ Date Last Updated:
+++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
+++ resultant, gcd
+++ Related Constructors: SparseMultivariatePolynomial, MultivariatePolynomial
+++ Also See:
+++ AMS Classifications:
+++ Keywords: polynomial, multivariate
+++ References:
+++ Description:
+++   This type is the basic representation of sparse recursive multivariate
+++ polynomials whose variables are arbitrary symbols. The ordering
+++ is alphabetic determined by the Symbol type.
+++ The coefficient ring may be non commutative,
+++ but the variables are assumed to commute.
+
+Polynomial(R:Ring):
+  PolynomialCategory(R, IndexedExponents Symbol, Symbol) with
+   if R has Algebra Fraction Integer then
+     integrate: (%, Symbol) -> %
+       ++ integrate(p,x) computes the integral of \spad{p*dx}, i.e.
+       ++ integrates the polynomial p with respect to the variable x.
+ == SparseMultivariatePolynomial(R, Symbol) add
+
+    import UserDefinedPartialOrdering(Symbol)
+
+    coerce(p:%):OutputForm ==
+      (r:= retractIfCan(p)@Union(R,"failed")) case R => r::R::OutputForm
+      a :=
+        userOrdered?() => largest variables p
+        mainVariable(p)::Symbol
+      outputForm(univariate(p, a), a::OutputForm)
+
+    if R has Algebra Fraction Integer then
+      integrate(p, x) == (integrate univariate(p, x)) (x::%)
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain IDEAL PolynomialIdeals}
 \pagehead{PolynomialIdeals}{IDEAL}
 \pagepic{ps/v103polynomialideals.ps}{IDEAL}{1.00}
@@ -39169,6 +43089,103 @@ RadicalFunctionField(F, UP, UPUP, radicnd, n): 
Exports == Impl where
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain RMATRIX RectangularMatrix}
+\pagehead{RectangularMatrix}{RMATRIX}
+\pagepic{ps/v103rectangularmatrix.ps}{RMATRIX}{1.00}
+See also:\\
+\refto{IndexedMatrix}{IMATRIX}
+\refto{Matrix}{MATRIX}
+\refto{SquareMatrix}{SQMATRIX}
+<<domain RMATRIX RectangularMatrix>>=
+)abbrev domain RMATRIX RectangularMatrix
+++ Author: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: IndexedMatrix, Matrix, SquareMatrix
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{RectangularMatrix} is a matrix domain where the number of rows
+++   and the number of columns are parameters of the domain.
+RectangularMatrix(m,n,R): Exports == Implementation where
+  m,n : NonNegativeInteger
+  R   : Ring
+  Row ==> DirectProduct(n,R)
+  Col ==> DirectProduct(m,R)
+  Exports ==> Join(RectangularMatrixCategory(m,n,R,Row,Col),_
+                   CoercibleTo Matrix R) with
+ 
+    if R has Field then VectorSpace R
+ 
+    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
+
+    rectangularMatrix: Matrix R -> $
+      ++ \spad{rectangularMatrix(m)} converts a matrix of type 
\spadtype{Matrix}
+      ++ to a matrix of type \spad{RectangularMatrix}.
+    coerce: $ -> Matrix R
+      ++ \spad{coerce(m)} converts a matrix of type 
\spadtype{RectangularMatrix}
+      ++ to a matrix of type \spad{Matrix}.
+ 
+  Implementation ==> Matrix R add
+    minr ==> minRowIndex
+    maxr ==> maxRowIndex
+    minc ==> minColIndex
+    maxc ==> maxColIndex
+    mini ==> minIndex
+    maxi ==> maxIndex
+ 
+    ZERO := new(m,n,0)$Matrix(R) pretend $
+    0    == ZERO
+ 
+    coerce(x:$):OutputForm == coerce(x pretend Matrix R)$Matrix(R)
+
+    matrix(l: List List R) ==
+      -- error check: this is a top level function
+      #l ^= m => error "matrix: wrong number of rows"
+      for ll in l repeat
+        #ll ^= n => error "matrix: wrong number of columns"
+      ans : Matrix R := new(m,n,0)
+      for i in minr(ans)..maxr(ans) for ll in l repeat
+        for j in minc(ans)..maxc(ans) for r in ll repeat
+          qsetelt_!(ans,i,j,r)
+      ans pretend $
+ 
+    row(x,i)    == directProduct row(x pretend Matrix(R),i)
+    column(x,j) == directProduct column(x pretend Matrix(R),j)
+ 
+    coerce(x:$):Matrix(R) == copy(x pretend Matrix(R))
+ 
+    rectangularMatrix x ==
+      (nrows(x) ^= m) or (ncols(x) ^= n) =>
+        error "rectangularMatrix: matrix of bad dimensions"
+      copy(x) pretend $
+ 
+    if R has EuclideanDomain then
+ 
+      rowEchelon x == rowEchelon(x pretend Matrix(R)) pretend $
+ 
+    if R has IntegralDomain then
+ 
+      rank x    == rank(x pretend Matrix(R))
+      nullity x == nullity(x pretend Matrix(R))
+      nullSpace x ==
+        [directProduct c for c in nullSpace(x pretend Matrix(R))]
+ 
+    if R has Field then
+ 
+      dimension() == (m * n) :: CardinalNumber
+ 
+    if R has ConvertibleTo InputForm then
+      convert(x:$):InputForm ==
+         convert [convert("rectangularMatrix"::Symbol)@InputForm,
+                  convert(x::Matrix(R))]$List(InputForm)
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain REF Reference}
 <<dot>>=
 "REF" -> "TYPE"
@@ -40584,6 +44601,878 @@ SimpleFortranProgram(R,FS): Exports == Implementation 
where
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain SAOS SingletonAsOrderedSet}
+<<dot>>=
+"SAOS" -> "ORDSET"
+"SingletonAsOrderedSet()" -> "OrderedSet()"
+@
+\pagehead{SingletonAsOrderedSet}{SAOS}
+\pagepic{ps/v103singletonasorderedset.ps}{SAOS}{1.00}
+<<domain SAOS SingletonAsOrderedSet>>=
+)abbrev domain SAOS SingletonAsOrderedSet
+++ This trivial domain lets us build Univariate Polynomials
+++  in an anonymous variable
+SingletonAsOrderedSet(): OrderedSet with
+              create:() -> %
+              convert:% -> Symbol
+  ==  add
+   create() == "?" pretend %
+   a<b == false -- only one element
+   coerce(a) == outputForm "?"  -- CJW doesn't like this: change ?
+   a=b == true  -- only one element
+   min(a,b) == a  -- only one element
+   max(a,b) == a  -- only one element
+   convert a == coerce("?")
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain SMP SparseMultivariatePolynomial}
+\pagehead{SparseMultivariatePolynomial}{SMP}
+\pagepic{ps/v103sparsemultivariatepolynomial.ps}{SMP}{1.00}
+See also:\\
+\refto{Polynomial}{POLY}
+\refto{MultivariatePolynomial}{MPOLY}
+\refto{IndexedExponents}{INDE}
+<<domain SMP SparseMultivariatePolynomial>>=
+)abbrev domain SMP SparseMultivariatePolynomial
+++ Author: Dave Barton, Barry Trager
+++ Date Created:
+++ Date Last Updated: 30 November 1994
+++ Fix History:
+++ 30 Nov 94: added gcdPolynomial for float-type coefficients
+++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
+++ resultant, gcd
+++ Related Constructors: Polynomial, MultivariatePolynomial
+++ Also See:
+++ AMS Classifications:
+++ Keywords: polynomial, multivariate
+++ References:
+++ Description:
+++   This type is the basic representation of sparse recursive multivariate
+++ polynomials. It is parameterized by the coefficient ring and the
+++ variable set which may be infinite. The variable ordering is determined
+++ by the variable set parameter. The coefficient ring may be non-commutative,
+++ but the variables are assumed to commute.
+
+SparseMultivariatePolynomial(R: Ring,VarSet: OrderedSet): C == T where
+  pgcd ==> PolynomialGcdPackage(IndexedExponents VarSet,VarSet,R,%)
+  C == PolynomialCategory(R,IndexedExponents(VarSet),VarSet)
+  SUP ==> SparseUnivariatePolynomial
+  T == add
+    --constants
+    --D := F(%) replaced by next line until compiler support completed
+
+    --representations
+      D := SparseUnivariatePolynomial(%)
+      VPoly:=  Record(v:VarSet,ts:D)
+      Rep:=  Union(R,VPoly)
+
+    --local function
+
+
+    --declarations
+      fn: R -> R
+      n: Integer
+      k: NonNegativeInteger
+      kp:PositiveInteger
+      k1:NonNegativeInteger
+      c: R
+      mvar: VarSet
+      val : R
+      var:VarSet
+      up: D
+      p,p1,p2,pval: %
+      Lval : List(R)
+      Lpval : List(%)
+      Lvar : List(VarSet)
+
+    --define
+      0  == 0$R::%
+      1  == 1$R::%
+
+
+      zero? p == p case R and zero?(p)$R
+--      one? p == p case R and one?(p)$R
+      one? p == p case R and ((p) = 1)$R
+    -- a local function
+      red(p:%):% ==
+         p case R => 0
+         if ground?(reductum p.ts) then leadingCoefficient(reductum p.ts) else 
[p.v,reductum p.ts]$VPoly
+
+      numberOfMonomials(p): NonNegativeInteger ==
+        p case R => 
+          zero?(p)$R => 0
+          1
+        +/[numberOfMonomials q for q in coefficients(p.ts)]
+
+      coerce(mvar):% == [mvar,monomial(1,1)$D]$VPoly
+
+      monomial? p ==
+        p case R => true
+        sup : D := p.ts
+        1 ^= numberOfMonomials(sup) => false
+        monomial? leadingCoefficient(sup)$D
+
+--    local
+      moreThanOneVariable?: % -> Boolean
+
+      moreThanOneVariable? p == 
+         p case R => false
+         q:=p.ts
+         any?(not ground? #1 ,coefficients q) => true
+         false
+
+      -- if we already know we use this (slighlty) faster function
+      univariateKnown: % -> SparseUnivariatePolynomial R 
+
+      univariateKnown p == 
+        p case R => (leadingCoefficient p) :: SparseUnivariatePolynomial(R)
+       monomial( leadingCoefficient p,degree p.ts)+ univariateKnown(red p)
+
+      univariate p ==
+        p case R =>(leadingCoefficient p) :: SparseUnivariatePolynomial(R)
+        moreThanOneVariable?  p => error "not univariate"
+        monomial( leadingCoefficient p,degree p.ts)+ univariate(red p)
+
+      multivariate (u:SparseUnivariatePolynomial(R),var:VarSet) ==
+        ground? u => (leadingCoefficient u) ::%
+        [var,monomial(leadingCoefficient u,degree u)$D]$VPoly +
+           multivariate(reductum u,var)
+
+      univariate(p:%,mvar:VarSet):SparseUnivariatePolynomial(%) ==
+        p case R or mvar>p.v  => monomial(p,0)$D
+        pt:=p.ts
+        mvar=p.v => pt
+        monomial(1,p.v,degree pt)*univariate(leadingCoefficient pt,mvar)+
+          univariate(red p,mvar)
+
+--  a local functions, used in next definition
+      unlikeUnivReconstruct(u:SparseUnivariatePolynomial(%),mvar:VarSet):% ==
+        zero? (d:=degree u) => coefficient(u,0)
+        monomial(leadingCoefficient u,mvar,d)+
+            unlikeUnivReconstruct(reductum u,mvar)
+
+      multivariate(u:SparseUnivariatePolynomial(%),mvar:VarSet):% ==
+        ground? u => coefficient(u,0)
+        uu:=u
+        while not zero? uu repeat
+          cc:=leadingCoefficient uu
+          cc case R or mvar > cc.v => uu:=reductum uu
+          return unlikeUnivReconstruct(u,mvar)
+        [mvar,u]$VPoly
+
+      ground?(p:%):Boolean ==
+        p case R => true
+        false
+
+--      const p ==
+--        p case R => p
+--        error "the polynomial is not a constant"
+
+      monomial(p,mvar,k1) ==
+        zero? k1 or zero? p => p
+        p case R or mvar>p.v => [mvar,monomial(p,k1)$D]$VPoly
+        p*[mvar,monomial(1,k1)$D]$VPoly
+
+      monomial(c:R,e:IndexedExponents(VarSet)):% ==
+        zero? e => (c::%)
+        monomial(1,leadingSupport e, leadingCoefficient e) *
+            monomial(c,reductum e)
+
+      coefficient(p:%, e:IndexedExponents(VarSet)) : R ==
+        zero? e =>
+          p case R  => p::R
+          coefficient(coefficient(p.ts,0),e)
+        p case R => 0
+        ve := leadingSupport e
+        vp := p.v
+        ve < vp =>
+          coefficient(coefficient(p.ts,0),e)
+        ve > vp => 0
+        coefficient(coefficient(p.ts,leadingCoefficient e),reductum e)
+
+--    coerce(e:IndexedExponents(VarSet)) : % ==
+--      e = 0 => 1
+--      monomial(1,leadingSupport e, leadingCoefficient e) *
+--          (reductum e)::%
+
+--    retract(p:%):IndexedExponents(VarSet) ==
+--      q:Union(IndexedExponents(VarSet),"failed"):=retractIfCan p
+--      q :: IndexedExponents(VarSet)
+
+--    retractIfCan(p:%):Union(IndexedExponents(VarSet),"failed") ==
+--      p = 0 => degree p
+--      reductum(p)=0 and leadingCoefficient(p)=1 => degree p
+--      "failed"
+
+      coerce(n) == n::R::%
+      coerce(c) == c::%
+      characteristic == characteristic$R
+
+      recip(p) ==
+        p case R => (uu:=recip(p::R);uu case "failed" => "failed"; uu::%)
+        "failed"
+
+      - p ==
+          p case R => -$R p
+          [p.v, - p.ts]$VPoly
+      n * p  ==
+          p case R => n * p::R
+          mvar:=p.v
+          up:=n*p.ts
+          if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+      c * p  ==
+          c = 1 => p
+          p case R => c * p::R
+          mvar:=p.v
+          up:=c*p.ts
+          if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+      p1 + p2  ==
+         p1 case R and p2 case R => p1 +$R p2
+         p1 case R => [p2.v, p1::D + p2.ts]$VPoly
+        p2 case R => [p1.v,  p1.ts + p2::D]$VPoly
+        p1.v = p2.v => 
+              mvar:=p1.v
+              up:=p1.ts+p2.ts
+              if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+        p1.v < p2.v =>
+              [p2.v, p1::D + p2.ts]$VPoly
+         [p1.v, p1.ts + p2::D]$VPoly
+
+      p1 - p2  ==
+         p1 case R and p2 case R => p1 -$R p2
+         p1 case R => [p2.v, p1::D - p2.ts]$VPoly
+         p2 case R => [p1.v,  p1.ts - p2::D]$VPoly
+         p1.v = p2.v =>
+              mvar:=p1.v
+              up:=p1.ts-p2.ts
+              if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+         p1.v < p2.v =>
+              [p2.v, p1::D - p2.ts]$VPoly
+         [p1.v, p1.ts - p2::D]$VPoly
+
+      p1 = p2  ==
+         p1 case R =>
+             p2 case R => p1 =$R p2
+             false
+         p2 case R => false
+         p1.v = p2.v => p1.ts = p2.ts
+         false
+
+      p1 * p2  ==
+         p1 case R => p1::R * p2
+         p2 case R => 
+            mvar:=p1.v
+            up:=p1.ts*p2
+            if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+        p1.v = p2.v => 
+            mvar:=p1.v
+            up:=p1.ts*p2.ts
+            if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+         p1.v > p2.v => 
+            mvar:=p1.v
+            up:=p1.ts*p2
+            if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+            --- p1.v < p2.v 
+         mvar:=p2.v
+         up:=p1*p2.ts
+         if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+
+      p ^ kp == p ** (kp pretend NonNegativeInteger)
+      p ** kp == p ** (kp pretend NonNegativeInteger )
+      p ^ k == p ** k
+      p ** k  ==
+         p case R => p::R ** k
+         -- univariate special case 
+         not moreThanOneVariable? p => 
+             multivariate( (univariateKnown p) ** k , p.v)
+         mvar:=p.v
+         up:=p.ts ** k
+         if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+
+      if R has IntegralDomain then
+         UnitCorrAssoc ==> Record(unit:%,canonical:%,associate:%)
+         unitNormal(p) ==
+            u,c,a:R
+            p case R =>
+              (u,c,a):= unitNormal(p::R)$R
+              [u::%,c::%,a::%]$UnitCorrAssoc
+            (u,c,a):= unitNormal(leadingCoefficient(p))$R
+            [u::%,(a*p)::%,a::%]$UnitCorrAssoc
+         unitCanonical(p) ==
+            p case R => unitCanonical(p::R)$R
+            (u,c,a):= unitNormal(leadingCoefficient(p))$R
+            a*p
+         unit? p ==
+            p case R => unit?(p::R)$R
+            false
+         associates?(p1,p2) ==
+            p1 case R => p2 case R and associates?(p1,p2)$R
+            p2 case VPoly and p1.v = p2.v and associates?(p1.ts,p2.ts)
+
+         if R has approximate then
+           p1  exquo  p2  ==
+              p1 case R and p2 case R =>
+                a:= (p1::R  exquo  p2::R)
+                if a case "failed" then "failed" else a::%
+              zero? p1 => p1
+--              one? p2 => p1
+              (p2 = 1) => p1
+              p1 case R or p2 case VPoly and p1.v < p2.v => "failed"
+              p2 case R or p1.v > p2.v =>
+                 a:= (p1.ts  exquo  p2::D)
+                 a case "failed" => "failed"
+                 [p1.v,a]$VPoly::%
+              -- The next test is useful in the case that R has inexact
+              -- arithmetic (in particular when it is Interval(...)).
+              -- In the case where the test succeeds, empirical evidence
+              -- suggests that it can speed up the computation several times,
+              -- but in other cases where there are a lot of variables
+              -- and p1 and p2 differ only in the low order terms (e.g. 
p1=p2+1)
+              -- it slows exquo down by about 15-20%.
+              p1 = p2 => 1
+              a:= p1.ts  exquo  p2.ts
+              a case "failed" => "failed"
+              mvar:=p1.v
+              up:SUP %:=a
+              if ground? (up) then leadingCoefficient(up) else 
[mvar,up]$VPoly::%
+         else
+           p1  exquo  p2  ==
+              p1 case R and p2 case R =>
+                a:= (p1::R  exquo  p2::R)
+                if a case "failed" then "failed" else a::%
+              zero? p1 => p1
+--              one? p2 => p1
+              (p2 = 1) => p1
+              p1 case R or p2 case VPoly and p1.v < p2.v => "failed"
+              p2 case R or p1.v > p2.v =>
+                 a:= (p1.ts  exquo  p2::D)
+                 a case "failed" => "failed"
+                 [p1.v,a]$VPoly::%
+              a:= p1.ts  exquo  p2.ts
+              a case "failed" => "failed"
+              mvar:=p1.v
+              up:SUP %:=a
+              if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly::%
+
+      map(fn,p) ==
+         p case R => fn(p)
+         mvar:=p.v
+         up:=map(map(fn,#1),p.ts)
+         if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+
+      if R has Field then
+        (p : %) / (r : R) == inv(r) * p
+
+      if R has GcdDomain then
+        content(p) ==
+           p case R => p
+           c :R :=0
+           up:=p.ts
+--           while not(zero? up) and not(one? c) repeat
+           while not(zero? up) and not(c = 1) repeat
+               c:=gcd(c,content leadingCoefficient(up))
+               up := reductum up
+           c
+
+      if R has EuclideanDomain and R has CharacteristicZero and not(R has 
FloatingPointSystem)  then
+        content(p,mvar) ==
+          p case R => p
+          gcd(coefficients univariate(p,mvar))$pgcd
+
+        gcd(p1,p2) ==  gcd(p1,p2)$pgcd
+
+        gcd(lp:List %) ==  gcd(lp)$pgcd
+
+        gcdPolynomial(a:SUP $,b:SUP $):SUP $ == gcd(a,b)$pgcd
+
+      else if R has GcdDomain then
+        content(p,mvar) ==
+          p case R => p
+          content univariate(p,mvar)
+
+        gcd(p1,p2) ==
+           p1 case R =>
+              p2 case R => gcd(p1,p2)$R::%
+              zero? p1 => p2
+              gcd(p1, content(p2.ts))
+           p2 case R =>
+              zero? p2 => p1
+              gcd(p2, content(p1.ts))
+           p1.v < p2.v => gcd(p1, content(p2.ts))
+           p1.v > p2.v => gcd(content(p1.ts), p2)
+           mvar:=p1.v
+           up:=gcd(p1.ts, p2.ts)
+           if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+
+        if R has FloatingPointSystem then
+           -- eventually need a better notion of gcd's over floats
+           -- this essentially computes the gcds of the monomial contents
+           gcdPolynomial(a:SUP $,b:SUP $):SUP $ ==
+             ground? (a) =>
+                  zero? a => b
+                 gcd(leadingCoefficient a, content b)::SUP $
+             ground?(b) =>
+                  zero? b => b
+                 gcd(leadingCoefficient b, content a)::SUP $
+             conta := content a
+             mona:SUP $ := monomial(conta, minimumDegree a)
+              if mona ^= 1 then
+                  a := (a exquo mona)::SUP $
+             contb := content b
+             monb:SUP $ := monomial(contb, minimumDegree b)
+              if monb ^= 1 then
+                  b := (b exquo monb)::SUP $
+             mong:SUP $  := monomial(gcd(conta, contb),
+                                      min(degree mona, degree monb))
+              degree(a) >= degree b =>
+                  not((a exquo b) case "failed") =>
+                       mong * b
+                  mong
+             not((b exquo a) case "failed") => mong * a
+             mong
+
+      coerce(p):OutputForm ==
+        p case R => (p::R)::OutputForm
+        outputForm(p.ts,p.v::OutputForm)
+
+      coefficients p ==
+        p case R => list(p :: R)$List(R)
+        "append"/[coefficients(p1)$% for p1 in coefficients(p.ts)]
+
+      retract(p:%):R ==
+        p case R => p :: R
+        error "cannot retract nonconstant polynomial"
+
+      retractIfCan(p:%):Union(R, "failed") ==
+        p case R => p::R
+        "failed"
+
+--      leadingCoefficientRecursive(p:%):% ==
+--         p case R => p
+--         leadingCoefficient p.ts
+
+      mymerge:(List VarSet,List VarSet) ->List VarSet
+      mymerge(l:List VarSet,m:List VarSet):List VarSet  ==
+         empty? l => m
+         empty? m => l
+         first l = first m => 
+            empty? rest l => 
+                 setrest!(l,rest m)
+                 l
+            empty? rest m => l 
+           setrest!(l, mymerge(rest l, rest m))
+            l
+         first l > first m =>
+            empty? rest l => 
+                setrest!(l,m) 
+                l
+            setrest!(l, mymerge(rest l, m))
+            l
+         empty? rest m => 
+             setrest!(m,l)
+             m
+         setrest!(m,mymerge(l,rest m))
+         m
+         
+      variables p ==
+         p case R => empty()
+         lv:List VarSet:=empty()
+         q := p.ts
+         while not zero? q repeat
+           lv:=mymerge(lv,variables leadingCoefficient q)
+           q := reductum q
+         cons(p.v,lv)
+
+      mainVariable p ==
+         p case R => "failed"
+         p.v
+
+      eval(p,mvar,pval) == univariate(p,mvar)(pval)
+      eval(p,mvar,val) ==  univariate(p,mvar)(val)
+
+      evalSortedVarlist(p,Lvar,Lpval):% ==
+        p case R => p
+        empty? Lvar or empty? Lpval => p
+        mvar := Lvar.first
+        mvar > p.v => evalSortedVarlist(p,Lvar.rest,Lpval.rest)
+        pval := Lpval.first
+        pts := map(evalSortedVarlist(#1,Lvar,Lpval),p.ts)
+        mvar=p.v =>
+             pval case R => pts (pval::R)
+             pts pval
+        multivariate(pts,p.v)
+
+      eval(p,Lvar,Lpval) ==
+       empty? rest Lvar => evalSortedVarlist(p,Lvar,Lpval)
+       sorted?(#1 > #2, Lvar) => evalSortedVarlist(p,Lvar,Lpval)
+        nlvar := sort(#1 > #2,Lvar)
+        nlpval :=
+           Lvar = nlvar => Lpval
+           nlpval := [Lpval.position(mvar,Lvar) for mvar in nlvar]
+        evalSortedVarlist(p,nlvar,nlpval)
+
+      eval(p,Lvar,Lval) ==
+        eval(p,Lvar,[val::% for val in Lval]$(List %)) -- kill?
+
+      degree(p,mvar) ==
+        p case R => 0
+        mvar= p.v => degree p.ts
+        mvar > p.v => 0    -- might as well take advantage of the order
+        max(degree(leadingCoefficient p.ts,mvar),degree(red p,mvar))
+
+      degree(p,Lvar)  == [degree(p,mvar)  for mvar in Lvar]
+
+      degree p ==
+        p case R => 0
+        degree(leadingCoefficient(p.ts)) + monomial(degree(p.ts), p.v)
+
+      minimumDegree p ==
+        p case R => 0
+        md := minimumDegree p.ts
+        minimumDegree(coefficient(p.ts,md)) + monomial(md, p.v)
+
+      minimumDegree(p,mvar) ==
+        p case R => 0
+        mvar = p.v => minimumDegree p.ts
+        md:=minimumDegree(leadingCoefficient p.ts,mvar)
+        zero? (p1:=red p) => md
+        min(md,minimumDegree(p1,mvar))
+
+      minimumDegree(p,Lvar) ==
+        [minimumDegree(p,mvar) for mvar in Lvar]
+
+      totalDegree(p, Lvar) ==
+        ground? p => 0
+        null setIntersection(Lvar, variables p) => 0
+        u := univariate(p, mv := mainVariable(p)::VarSet)
+        weight:NonNegativeInteger := (member?(mv,Lvar) => 1; 0)
+        tdeg:NonNegativeInteger := 0
+        while u ^= 0 repeat
+            termdeg:NonNegativeInteger := weight*degree u +
+                           totalDegree(leadingCoefficient u, Lvar)
+            tdeg := max(tdeg, termdeg)
+            u := reductum u
+        tdeg
+
+      if R has CommutativeRing then
+        differentiate(p,mvar) ==
+          p case R => 0
+          mvar=p.v =>  
+             up:=differentiate p.ts
+             if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
+          up:=map(differentiate(#1,mvar),p.ts)
+          if ground? up then leadingCoefficient(up) else [p.v,up]$VPoly
+
+      leadingCoefficient(p) ==
+         p case R => p
+         leadingCoefficient(leadingCoefficient(p.ts))
+
+--      trailingCoef(p) ==
+--        p case R => p
+--        coef(p.ts,0) case R => coef(p.ts,0)
+--        trailingCoef(red p)
+--      TrailingCoef(p) == trailingCoef(p)
+
+      leadingMonomial p ==
+          p case R => p
+          monomial(leadingMonomial leadingCoefficient(p.ts),
+                   p.v, degree(p.ts))
+
+      reductum(p) == 
+          p case R => 0
+          p - leadingMonomial p
+
+
+--        if R is Integer then
+--           pgcd := PolynomialGcdPackage(%,VarSet)
+--           gcd(p1,p2) ==
+--               gcd(p1,p2)$pgcd
+--
+--        else if R is RationalNumber then
+--           gcd(p1,p2) ==
+--               mrat:= MRationalFactorize(VarSet,%)
+--               gcd(p1,p2)$mrat
+--
+--        else gcd(p1,p2) ==
+--           p1 case R =>
+--              p2 case R => gcd(p1,p2)$R::%
+--              p1 = 0 => p2
+--              gcd(p1, content(p2.ts))
+--           p2 case R =>
+--              p2 = 0 => p1
+--              gcd(p2, content(p1.ts))
+--           p1.v < p2.v => gcd(p1, content(p2.ts))
+--           p1.v > p2.v => gcd(content(p1.ts), p2)
+--           PSimp(p1.v, gcd(p1.ts, p2.ts))
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain SMTS SparseMultivariateTaylorSeries}
+\pagehead{SparseMultivariateTaylorSeries}{SMTS}
+\pagepic{ps/v103sparsemultivariatetaylorseries.ps}{SMTS}{1.00}
+See also:\\
+\refto{TaylorSeries}{TS}
+<<domain SMTS SparseMultivariateTaylorSeries>>=
+)abbrev domain SMTS SparseMultivariateTaylorSeries
+++ This domain provides multivariate Taylor series
+++ Authors: William Burge, Stephen Watt, Clifton Williamson
+++ Date Created: 15 August 1988
+++ Date Last Updated: 18 May 1991
+++ Basic Operations:
+++ Related Domains:
+++ Also See: UnivariateTaylorSeries
+++ AMS Classifications:
+++ Keywords: multivariate, Taylor, series
+++ Examples:
+++ References:
+++ Description:
+++   This domain provides multivariate Taylor series with variables
+++   from an arbitrary ordered set.  A Taylor series is represented
+++   by a stream of polynomials from the polynomial domain SMP.
+++   The nth element of the stream is a form of degree n.  SMTS is an
+++   internal domain.
+SparseMultivariateTaylorSeries(Coef,Var,SMP):_
+ Exports == Implementation where
+  Coef : Ring
+  Var  : OrderedSet
+  SMP  : PolynomialCategory(Coef,IndexedExponents Var,Var)
+  I   ==> Integer
+  L   ==> List
+  NNI ==> NonNegativeInteger
+  OUT ==> OutputForm
+  PS  ==> InnerTaylorSeries SMP
+  RN  ==> Fraction Integer
+  ST  ==> Stream
+  StS ==> Stream SMP
+  STT ==> StreamTaylorSeriesOperations SMP
+  STF ==> StreamTranscendentalFunctions SMP
+  ST2 ==> StreamFunctions2
+  ST3 ==> StreamFunctions3
+ 
+  Exports ==> MultivariateTaylorSeriesCategory(Coef,Var) with
+    coefficient: (%,NNI) -> SMP
+      ++ \spad{coefficient(s, n)} gives the terms of total degree n.
+    coerce: Var -> %
+      ++ \spad{coerce(var)} converts a variable to a Taylor series
+    coerce: SMP -> %
+      ++ \spad{coerce(poly)} regroups the terms by total degree and forms
+      ++ a series.
+    "*":(SMP,%)->%
+      ++\spad{smp*ts} multiplies a TaylorSeries by a monomial SMP.
+    csubst:(L Var,L StS) -> (SMP -> StS)
+      ++\spad{csubst(a,b)} is for internal use only
+ 
+    if Coef has Algebra Fraction Integer then
+      integrate: (%,Var,Coef) -> %
+        ++\spad{integrate(s,v,c)} is the integral of s with respect
+        ++ to v and having c as the constant of integration.
+      fintegrate: (() -> %,Var,Coef) -> %
+        ++\spad{fintegrate(f,v,c)} is the integral of \spad{f()} with respect
+        ++ to v and having c as the constant of integration.
+        ++ The evaluation of \spad{f()} is delayed.
+ 
+  Implementation ==> PS add
+ 
+    Rep := StS -- Below we use the fact that Rep of PS is Stream SMP.
+    extend(x,n) == extend(x,n + 1)$Rep
+    complete x == complete(x)$Rep
+
+    evalstream:(%,L Var,L SMP) -> StS
+    evalstream(s:%,lv:(L Var),lsmp:(L SMP)):(ST SMP) ==
+      scan(0,_+$SMP,map(eval(#1,lv,lsmp),s pretend StS))$ST2(SMP,SMP)
+ 
+    addvariable:(Var,InnerTaylorSeries Coef) -> %
+    addvariable(v,s) ==
+      ints := integers(0)$STT pretend ST NNI
+      map(monomial(#2 :: SMP,v,#1)$SMP,ints,s pretend ST 
Coef)$ST3(NNI,Coef,SMP)
+ 
+    coefficient(s,n) == elt(s,n + 1)$Rep  -- 1-based indexing for streams
+ 
+--% creation of series
+ 
+    coerce(r:Coef) == monom(r::SMP,0)$STT
+    smp:SMP * p:% == (((smp) *  (p pretend Rep))$STT)pretend %
+    r:Coef * p:% == (((r::SMP) *  (p pretend Rep))$STT)pretend %
+    p:% * r:Coef == (((r::SMP) * ( p pretend Rep))$STT)pretend %
+    mts(p:SMP):% ==
+      (uv := mainVariable p) case "failed" => monom(p,0)$STT
+      v := uv :: Var
+      s : % := 0
+      up := univariate(p,v)
+      while not zero? up repeat
+        s := s + monomial(1,v,degree up) * mts(leadingCoefficient up)
+        up := reductum up
+      s
+ 
+    coerce(p:SMP) == mts p
+    coerce(v:Var) == v :: SMP :: %
+ 
+    monomial(r:%,v:Var,n:NNI) ==
+      r * monom(monomial(1,v,n)$SMP,n)$STT
+ 
+--% evaluation
+ 
+    substvar: (SMP,L Var,L %) -> %
+    substvar(p,vl,q) ==
+      null vl => monom(p,0)$STT
+      (uv := mainVariable p) case "failed" => monom(p,0)$STT
+      v := uv :: Var
+      v = first vl =>
+        s : % := 0
+        up := univariate(p,v)
+        while not zero? up repeat
+          c := leadingCoefficient up
+          s := s + first q ** degree up * substvar(c,rest vl,rest q)
+          up := reductum up
+        s
+      substvar(p,rest vl,rest q)
+ 
+    sortmfirst:(SMP,L Var,L %) -> %
+    sortmfirst(p,vl,q) ==
+      nlv : L Var := sort(#1 > #2,vl)
+      nq : L % := [q position$(L Var) (i,vl) for i in nlv]
+      substvar(p,nlv,nq)
+ 
+    csubst(vl,q) == sortmfirst(#1,vl,q pretend L(%)) pretend StS
+ 
+    restCheck(s:StS):StS ==
+      -- checks that stream is null or first element is 0
+      -- returns empty() or rest of stream
+      empty? s => s
+      not zero? frst s =>
+        error "eval: constant coefficient should be 0"
+      rst s
+ 
+    eval(s:%,v:L Var,q:L %) ==
+      #v ^= #q =>
+        error "eval: number of variables should equal number of values"
+      nq : L StS := [restCheck(i pretend StS) for i in q]
+      addiag(map(csubst(v,nq),s pretend StS)$ST2(SMP,StS))$STT pretend %
+ 
+    substmts(v:Var,p:SMP,q:%):% ==
+      up := univariate(p,v)
+      ss : % := 0
+      while not zero? up repeat
+        d:=degree up
+        c:SMP:=leadingCoefficient up
+        ss := ss + c* q ** d
+        up := reductum up
+      ss
+ 
+    subststream(v:Var,p:SMP,q:StS):StS==
+      substmts(v,p,q pretend %) pretend StS
+ 
+    comp1:(Var,StS,StS) -> StS
+    comp1(v,r,t)== addiag(map(subststream(v,#1,t),r)$ST2(SMP,StS))$STT
+ 
+    comp(v:Var,s:StS,t:StS):StS == delay
+      empty? s => s
+      f := frst s; r : StS := rst s;
+      empty? r => s
+      empty? t => concat(f,comp1(v,r,empty()$StS))
+      not zero? frst t =>
+        error "eval: constant coefficient should be zero"
+      concat(f,comp1(v,r,rst t))
+ 
+    eval(s:%,v:Var,t:%) == comp(v,s pretend StS,t pretend StS)
+ 
+--% differentiation and integration
+ 
+    differentiate(s:%,v:Var):% ==
+      empty? s => 0
+      map(differentiate(#1,v),rst s)
+ 
+    if Coef has Algebra Fraction Integer then
+ 
+      stream(x:%):Rep == x pretend Rep
+ 
+      (x:%) ** (r:RN) == powern(r,stream x)$STT
+      (r:RN) * (x:%)  == map(r * #1, stream x)$ST2(SMP,SMP) pretend %
+      (x:%) * (r:RN)  == map(#1 * r,stream x )$ST2(SMP,SMP) pretend %
+ 
+      exp x == exp(stream x)$STF
+      log x == log(stream x)$STF
+ 
+      sin x == sin(stream x)$STF
+      cos x == cos(stream x)$STF
+      tan x == tan(stream x)$STF
+      cot x == cot(stream x)$STF
+      sec x == sec(stream x)$STF
+      csc x == csc(stream x)$STF
+ 
+      asin x == asin(stream x)$STF
+      acos x == acos(stream x)$STF
+      atan x == atan(stream x)$STF
+      acot x == acot(stream x)$STF
+      asec x == asec(stream x)$STF
+      acsc x == acsc(stream x)$STF
+ 
+      sinh x == sinh(stream x)$STF
+      cosh x == cosh(stream x)$STF
+      tanh x == tanh(stream x)$STF
+      coth x == coth(stream x)$STF
+      sech x == sech(stream x)$STF
+      csch x == csch(stream x)$STF
+ 
+      asinh x == asinh(stream x)$STF
+      acosh x == acosh(stream x)$STF
+      atanh x == atanh(stream x)$STF
+      acoth x == acoth(stream x)$STF
+      asech x == asech(stream x)$STF
+      acsch x == acsch(stream x)$STF
+ 
+      intsmp(v:Var,p: SMP): SMP ==
+        up := univariate(p,v)
+        ss : SMP := 0
+        while not zero? up repeat
+          d := degree up
+          c := leadingCoefficient up
+          ss := ss + inv((d+1) :: RN) * monomial(c,v,d+1)$SMP
+          up := reductum up
+        ss
+ 
+      fintegrate(f,v,r) ==
+        concat(r::SMP,delay map(intsmp(v,#1),f() pretend StS))
+      integrate(s,v,r) ==
+        concat(r::SMP,map(intsmp(v,#1),s pretend StS))
+ 
+    -- If there is more than one term of the same order, group them.
+    tout(p:SMP):OUT ==
+      pe := p :: OUT
+      monomial? p => pe
+      paren pe
+ 
+    showAll?: () -> Boolean
+    -- check a global Lisp variable
+    showAll?() == true
+ 
+    coerce(s:%):OUT ==
+      uu := s pretend Stream(SMP)
+      empty? uu => (0$SMP) :: OUT
+      n : NNI; count : NNI := _$streamCount$Lisp
+      l : List OUT := empty()
+      for n in 0..count while not empty? uu repeat
+        if frst(uu) ^= 0 then l := concat(tout frst uu,l)
+        uu := rst uu
+      if showAll?() then
+        for n in n.. while explicitEntries? uu and _
+               not eq?(uu,rst uu) repeat
+          if frst(uu) ^= 0 then l := concat(tout frst uu,l)
+          uu := rst uu
+      l :=
+        explicitlyEmpty? uu => l
+        eq?(uu,rst uu) and frst uu = 0 => l
+        concat(prefix("O" :: OUT,[n :: OUT]),l)
+      empty? l => (0$SMP) :: OUT
+      reduce("+",reverse_! l)
+    if Coef has Field then
+         stream(x:%):Rep == x pretend Rep
+         SF2==> StreamFunctions2
+         p:% / r:Coef ==(map(#1/$SMP r,stream p)$SF2(SMP,SMP))pretend %
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain SHDP SplitHomogeneousDirectProduct}
 \pagehead{SplitHomogeneousDirectProduct}{SHDP}
 \pagepic{ps/v103splithomogeneousdirectproduct.ps}{SHDP}{1.00}
@@ -40640,6 +45529,296 @@ SplitHomogeneousDirectProduct(dimtot,dim1,S) : T == C 
where
 
 @
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain SQMATRIX SquareMatrix}
+<<SquareMatrix.input>>=
+-- matrix.spad.pamphlet SquareMatrix.input
+)spool SquareMatrix.output
+)set message test on
+)set message auto off
+)clear all
+--S 1 of 6
+)set expose add constructor SquareMatrix
+--R 
+--I   SquareMatrix is now explicitly exposed in frame frame0 
+--E 1 
+
+--S 2 of 6
+m := squareMatrix [ [1,-%i],[%i,4] ]
+--R 
+--R
+--R        +1   - %i+
+--R   (1)  |        |
+--R        +%i   4  +
+--R                                        Type: SquareMatrix(2,Complex 
Integer)
+--E 2
+
+--S 3 of 6
+m*m - m
+--R 
+--R
+--R        + 1   - 4%i+
+--R   (2)  |          |
+--R        +4%i   13  +
+--R                                        Type: SquareMatrix(2,Complex 
Integer)
+--E 3
+
+--S 4 of 6
+mm := squareMatrix [ [m, 1], [1-m, m**2] ]
+--R 
+--R
+--R        ++1   - %i+      +1  0+   +
+--R        ||        |      |    |   |
+--R        |+%i   4  +      +0  1+   |
+--R   (3)  |                         |
+--R        |+ 0    %i +  + 2   - 5%i+|
+--R        ||         |  |          ||
+--R        ++- %i  - 3+  +5%i   17  ++
+--R                        Type: SquareMatrix(2,SquareMatrix(2,Complex 
Integer))
+--E 4
+
+--S 5 of 6
+p := (x + m)**2
+--R 
+--R
+--R         2   + 2   - 2%i+    + 2   - 5%i+
+--R   (4)  x  + |          |x + |          |
+--R             +2%i    8  +    +5%i   17  +
+--R                             Type: Polynomial SquareMatrix(2,Complex 
Integer)
+--E 5
+
+--S 6 of 6
+p::SquareMatrix(2, ?)
+--R 
+--R
+--R        + 2                        +
+--R        |x  + 2x + 2  - 2%i x - 5%i|
+--R   (5)  |                          |
+--R        |              2           |
+--R        +2%i x + 5%i  x  + 8x + 17 +
+--R                             Type: SquareMatrix(2,Polynomial Complex 
Integer)
+--E 6
+)spool
+)lisp (bye)
+@
+<<SquareMatrix.help>>=
+====================================================================
+SquareMatrix examples
+====================================================================
+ 
+The top level matrix type in Axiom is Matrix, which provides basic
+arithmetic and linear algebra functions.  However, since the matrices
+can be of any size it is not true that any pair can be added or
+multiplied.  Thus Matrix has little algebraic structure.
+ 
+Sometimes you want to use matrices as coefficients for polynomials or
+in other algebraic contexts.  In this case, SquareMatrix should be
+used.  The domain SquareMatrix(n,R) gives the ring of n by n square
+matrices over R.
+ 
+Since SquareMatrix is not normally exposed at the top level, you must
+expose it before it can be used.
+
+  )set expose add constructor SquareMatrix
+
+Once SQMATRIX has been exposed, values can be created using the
+squareMatrix function.
+
+  m := squareMatrix [ [1,-%i],[%i,4] ]
+    +1   - %i+
+    |        |
+    +%i   4  +
+                        Type: SquareMatrix(2,Complex Integer)
+
+The usual arithmetic operations are available.
+
+  m*m - m
+    + 1   - 4%i+
+    |          |
+    +4%i   13  +
+                        Type: SquareMatrix(2,Complex Integer)
+
+Square matrices can be used where ring elements are required.
+For example, here is a matrix with matrix entries.
+
+  mm := squareMatrix [ [m, 1], [1-m, m**2] ]
+    ++1   - %i+      +1  0+   +
+    ||        |      |    |   |
+    |+%i   4  +      +0  1+   |
+    |                         |
+    |+ 0    %i +  + 2   - 5%i+|
+    ||         |  |          ||
+    ++- %i  - 3+  +5%i   17  ++
+                        Type: SquareMatrix(2,SquareMatrix(2,Complex Integer))
+
+Or you can construct a polynomial with  square matrix coefficients.
+
+  p := (x + m)**2
+     2   + 2   - 2%i+    + 2   - 5%i+
+    x  + |          |x + |          |
+         +2%i    8  +    +5%i   17  +
+                        Type: Polynomial SquareMatrix(2,Complex Integer)
+
+This value can be converted to a square matrix with polynomial coefficients.
+
+  p::SquareMatrix(2, ?)
+    + 2                        +
+    |x  + 2x + 2  - 2%i x - 5%i|
+    |                          |
+    |              2           |
+    +2%i x + 5%i  x  + 8x + 17 +
+                        Type: SquareMatrix(2,Polynomial Complex Integer)
+ 
+See Also:
+o )help Matrix
+o )show SquareMatrix
+o $AXIOM/doc/src/algebra/matrix.spad.dvi
+
+@
+\pagehead{SquareMatrix}{SQMATRIX}
+\pagepic{ps/v103squarematrix.ps}{SQMATRIX}{1.00}
+See also:\\
+\refto{IndexedMatrix}{IMATRIX}
+\refto{Matrix}{MATRIX}
+\refto{RectangularMatrix}{RMATRIX}
+<<domain SQMATRIX SquareMatrix>>=
+)abbrev domain SQMATRIX SquareMatrix
+++ Author: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: IndexedMatrix, Matrix, RectangularMatrix
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{SquareMatrix} is a matrix domain of square matrices, where the
+++   number of rows (= number of columns) is a parameter of the type.
+SquareMatrix(ndim,R): Exports == Implementation where
+  ndim : NonNegativeInteger
+  R    : Ring
+  Row ==> DirectProduct(ndim,R)
+  Col ==> DirectProduct(ndim,R)
+  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,$)
+ 
+  Exports ==> Join(SquareMatrixCategory(ndim,R,Row,Col),_
+                   CoercibleTo Matrix R) with
+ 
+    transpose: $ -> $
+      ++ \spad{transpose(m)} returns the transpose of the matrix m.
+    squareMatrix: Matrix R -> $
+      ++ \spad{squareMatrix(m)} converts a matrix of type \spadtype{Matrix}
+      ++ to a matrix of type \spadtype{SquareMatrix}.
+    coerce: $ -> Matrix R
+      ++ \spad{coerce(m)} converts a matrix of type \spadtype{SquareMatrix}
+      ++ to a matrix of type \spadtype{Matrix}.
+--  symdecomp : $ -> Record(sym:$,antisym:$)
+--    ++ \spad{symdecomp(m)} decomposes the matrix m as a sum of a symmetric
+--    ++ matrix \spad{m1} and an antisymmetric matrix \spad{m2}. The object
+--    ++ returned is the Record \spad{[m1,m2]}
+--  if R has commutative("*") then
+--    minorsVect: -> Vector(Union(R,"uncomputed")) --range: 1..2**n-1
+--      ++ \spad{minorsVect(m)} returns a vector of the minors of the matrix m
+    if R has commutative("*") then central
+      ++ the elements of the Ring R, viewed as diagonal matrices, commute
+      ++ with all matrices and, indeed, are the only matrices which commute
+      ++ with all matrices.
+    if R has commutative("*") and R has unitsKnown then unitsKnown
+      ++ the invertible matrices are simply the matrices whose determinants
+      ++ are units in the Ring R.
+    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
+ 
+  Implementation ==> Matrix R add
+    minr ==> minRowIndex
+    maxr ==> maxRowIndex
+    minc ==> minColIndex
+    maxc ==> maxColIndex
+    mini ==> minIndex
+    maxi ==> maxIndex
+ 
+    ZERO := scalarMatrix 0
+    0    == ZERO
+    ONE  := scalarMatrix 1
+    1    == ONE
+
+    characteristic() == characteristic()$R
+ 
+    matrix(l: List List R) ==
+      -- error check: this is a top level function
+      #l ^= ndim => error "matrix: wrong number of rows"
+      for ll in l repeat
+        #ll ^= ndim => error "matrix: wrong number of columns"
+      ans : Matrix R := new(ndim,ndim,0)
+      for i in minr(ans)..maxr(ans) for ll in l repeat
+        for j in minc(ans)..maxc(ans) for r in ll repeat
+          qsetelt_!(ans,i,j,r)
+      ans pretend $
+ 
+    row(x,i)    == directProduct row(x pretend Matrix(R),i)
+    column(x,j) == directProduct column(x pretend Matrix(R),j)
+    coerce(x:$):OutputForm == coerce(x pretend Matrix R)$Matrix(R)
+ 
+    scalarMatrix r == scalarMatrix(ndim,r)$Matrix(R) pretend $
+ 
+    diagonalMatrix l ==
+      #l ^= ndim =>
+        error "diagonalMatrix: wrong number of entries in list"
+      diagonalMatrix(l)$Matrix(R) pretend $
+ 
+    coerce(x:$):Matrix(R) == copy(x pretend Matrix(R))
+ 
+    squareMatrix x ==
+      (nrows(x) ^= ndim) or (ncols(x) ^= ndim) =>
+        error "squareMatrix: matrix of bad dimensions"
+      copy(x) pretend $
+ 
+    x:$ * v:Col ==
+      directProduct((x pretend Matrix(R)) * (v :: Vector(R)))
+ 
+    v:Row * x:$ ==
+      directProduct((v :: Vector(R)) * (x pretend Matrix(R)))
+ 
+    x:$ ** n:NonNegativeInteger ==
+      ((x pretend Matrix(R)) ** n) pretend $
+ 
+    if R has commutative("*") then
+ 
+      determinant x == determinant(x pretend Matrix(R))
+      minordet x    == minordet(x pretend Matrix(R))
+ 
+    if R has EuclideanDomain then
+ 
+      rowEchelon x == rowEchelon(x pretend Matrix(R)) pretend $
+ 
+    if R has IntegralDomain then
+ 
+      rank x    == rank(x pretend Matrix(R))
+      nullity x == nullity(x pretend Matrix(R))
+      nullSpace x ==
+        [directProduct c for c in nullSpace(x pretend Matrix(R))]
+ 
+    if R has Field then
+ 
+      dimension() == (m * n) :: CardinalNumber
+ 
+      inverse x ==
+        (u := inverse(x pretend Matrix(R))) case "failed" => "failed"
+        (u :: Matrix(R)) pretend $
+ 
+      x:$ ** n:Integer ==
+        ((x pretend Matrix(R)) ** n) pretend $
+ 
+      recip x == inverse x
+ 
+    if R has ConvertibleTo InputForm then
+      convert(x:$):InputForm ==
+         convert [convert("squareMatrix"::Symbol)@InputForm,
+                  convert(x::Matrix(R))]$List(InputForm)
+
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain STACK Stack}
 <<dot>>=
 "STACK" -> "SKAGG"
@@ -40993,6 +46172,63 @@ SymbolTable() : exports == implementation where
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \chapter{Chapter T}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain TS TaylorSeries}
+\pagehead{TaylorSeries}{TS}
+\pagepic{ps/v103taylorseries.ps}{TS}{1.00}
+See also:\\
+\refto{SparseMultivariateTaylorSeries}{SMTS}
+<<domain TS TaylorSeries>>=
+)abbrev domain TS TaylorSeries
+++ Authors: Burge, Watt, Williamson
+++ Date Created: 15 August 1988
+++ Date Last Updated: 18 May 1991
+++ Basic Operations:
+++ Related Domains: SparseMultivariateTaylorSeries
+++ Also See: UnivariateTaylorSeries
+++ AMS Classifications:
+++ Keywords: multivariate, Taylor, series
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{TaylorSeries} is a general multivariate Taylor series domain
+++   over the ring Coef and with variables of type Symbol.
+TaylorSeries(Coef): Exports == Implementation where
+  Coef  : Ring
+  L   ==> List
+  NNI ==> NonNegativeInteger
+  SMP ==> Polynomial Coef
+  StS ==> Stream SMP
+ 
+  Exports ==> MultivariateTaylorSeriesCategory(Coef,Symbol) with
+    coefficient: (%,NNI) -> SMP
+      ++\spad{coefficient(s, n)} gives the terms of total degree n.
+    coerce: Symbol -> %
+      ++\spad{coerce(s)} converts a variable to a Taylor series
+    coerce: SMP -> %
+      ++\spad{coerce(s)} regroups terms of s by total degree
+      ++ and forms a series.
+ 
+    if Coef has Algebra Fraction Integer then
+      integrate: (%,Symbol,Coef) -> %
+        ++\spad{integrate(s,v,c)} is the integral of s with respect
+        ++ to v and having c as the constant of integration.
+      fintegrate: (() -> %,Symbol,Coef) -> %
+        ++\spad{fintegrate(f,v,c)} is the integral of \spad{f()} with respect
+        ++ to v and having c as the constant of integration.
+        ++ The evaluation of \spad{f()} is delayed.
+ 
+  Implementation ==> SparseMultivariateTaylorSeries(Coef,Symbol,SMP) add
+    Rep := StS -- Below we use the fact that Rep of PS is Stream SMP.
+ 
+    polynomial(s,n) ==
+      sum : SMP := 0
+      for i in 0..n while not empty? s repeat
+        sum := sum + frst s
+        s:= rst s
+      sum
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{domain TEXTFILE TextFile}
 <<TextFile.input>>=
 -- files.spad.pamphlet TextFile.input
@@ -45407,6 +50643,7 @@ Note that this code is not included in the generated 
catdef.spad file.
 <<domain EQ Equation>>
 <<domain EXPEXPAN ExponentialExpansion>>
 <<domain EXPUPXS ExponentialOfUnivariatePuiseuxSeries>>
+<<domain EMR EuclideanModularRing>>
 <<domain EXPR Expression>>
 <<domain E04DGFA e04dgfAnnaType>>
 <<domain E04FDFA e04fdfAnnaType>>
@@ -45450,6 +50687,7 @@ Note that this code is not included in the generated 
catdef.spad file.
 <<domain FPARFRAC FullPartialFractionExpansion>>
 
 <<domain GDMP GeneralDistributedMultivariatePolynomial>>
+<<domain GMODPOL GeneralModulePolynomial>>
 <<domain GCNAALG GenericNonAssociativeAlgebra>>
 <<domain GSERIES GeneralUnivariatePowerSeries>>
 
@@ -45465,8 +50703,10 @@ Note that this code is not included in the generated 
catdef.spad file.
 <<domain IDPO IndexedDirectProductObject>>
 <<domain IDPOAM IndexedDirectProductOrderedAbelianMonoid>>
 <<domain IDPOAMS IndexedDirectProductOrderedAbelianMonoidSup>>
+<<domain INDE IndexedExponents>>
 <<domain IFARRAY IndexedFlexibleArray>>
 <<domain ILIST IndexedList>>
+<<domain IMATRIX IndexedMatrix>>
 <<domain IARRAY1 IndexedOneDimensionalArray>>
 <<domain IARRAY2 IndexedTwoDimensionalArray>>
 <<domain ITUPLE InfiniteTuple>>
@@ -45474,6 +50714,7 @@ Note that this code is not included in the generated 
catdef.spad file.
 <<domain IFF InnerFiniteField>>
 <<domain IFAMON InnerFreeAbelianMonoid>>
 <<domain IIARRAY2 InnerIndexedTwoDimensionalArray>>
+<<domain INFORM InputForm>>
 <<domain INT Integer>>
 <<domain ZMOD IntegerMod>>
 <<domain INTFTBL IntegrationFunctionsTable>>
@@ -45499,6 +50740,15 @@ Note that this code is not included in the generated 
catdef.spad file.
 <<domain MFLOAT MachineFloat>>
 <<domain MINT MachineInteger>>
 <<domain MKCHSET MakeCachableSet>>
+<<domain MATRIX Matrix>>
+<<domain MODMON ModMonic>>
+<<domain MODMONOM ModuleMonomial>>
+<<domain MODFIELD ModularField>>
+<<domain MODRING ModularRing>>
+<<domain MOEBIUS MoebiusTransform>>
+<<domain MRING MonoidRing>>
+<<domain MSET Multiset>>
+<<domain MPOLY MultivariatePolynomial>>
 
 <<domain NONE None>>
 <<domain NNI NonNegativeInteger>>
@@ -45521,6 +50771,7 @@ Note that this code is not included in the generated 
catdef.spad file.
 <<domain ACPLOT PlaneAlgebraicCurvePlot>>
 <<domain PALETTE Palette>>
 <<domain HACKPI Pi>>
+<<domain POLY Polynomial>>
 <<domain IDEAL PolynomialIdeals>>
 <<domain PI PositiveInteger>>
 <<domain PRIMARR PrimitiveArray>>
@@ -45530,6 +50781,7 @@ Note that this code is not included in the generated 
catdef.spad file.
 <<domain QUEUE Queue>>
 
 <<domain RADFF RadicalFunctionField>>
+<<domain RMATRIX RectangularMatrix>>
 <<domain REF Reference>>
 <<domain RESULT Result>>
 <<domain ROMAN RomanNumeral>>
@@ -45537,14 +50789,19 @@ Note that this code is not included in the generated 
catdef.spad file.
 <<domain FORMULA ScriptFormulaFormat>>
 <<domain SAE SimpleAlgebraicExtension>>
 <<domain SFORT SimpleFortranProgram>>
+<<domain SAOS SingletonAsOrderedSet>>
 <<domain SDPOL SequentialDifferentialPolynomial>>
 <<domain SDVAR SequentialDifferentialVariable>>
 <<domain SETMN SetOfMIntegersInOneToN>>
+<<domain SMP SparseMultivariatePolynomial>>
+<<domain SMTS SparseMultivariateTaylorSeries>>
 <<domain SHDP SplitHomogeneousDirectProduct>>
+<<domain SQMATRIX SquareMatrix>>
 <<domain STACK Stack>>
 <<domain SWITCH Switch>>
 <<domain SYMTAB SymbolTable>>
 
+<<domain TS TaylorSeries>>
 <<domain TEXTFILE TextFile>>
 <<domain SYMS TheSymbolTable>>
 <<domain M3D ThreeDimensionalMatrix>>
diff --git a/src/algebra/acplot.spad.pamphlet b/src/algebra/acplot.spad.pamphlet
index 4dbec57..f18549a 100644
--- a/src/algebra/acplot.spad.pamphlet
+++ b/src/algebra/acplot.spad.pamphlet
@@ -125,9 +125,7 @@ ans3 := realSolve(lp,lsv,0.01)$REALSOLV
 --R                                                        Type: List List 
Float
 --E 13
 )spool
- 
-)spool
-)lisp (bye)
+
 @
 <<RealSolvePackage.help>>=
 ====================================================================
diff --git a/src/algebra/mappkg.spad.pamphlet b/src/algebra/mappkg.spad.pamphlet
index 030d679..ce3616e 100644
--- a/src/algebra/mappkg.spad.pamphlet
+++ b/src/algebra/mappkg.spad.pamphlet
@@ -318,14 +318,6 @@ fibs := curry(shiftfib, fibinit)
 )lisp (bye)
  
 @
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
-)spool
-)lisp (bye)
-@
 <<MappingPackage1.help>>=
 ====================================================================
 MappingPackage examples
@@ -1690,14 +1682,6 @@ h:=(x:EXPR(INT)):EXPR(INT)+->1
 )lisp (bye)
  
 @
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
-)spool
-)lisp (bye)
-@
 <<MappingPackage4.help>>=
 ====================================================================
 MappingPackage examples
diff --git a/src/algebra/xlpoly.spad.pamphlet b/src/algebra/xlpoly.spad.pamphlet
index 43702dd..bc703ed 100644
--- a/src/algebra/xlpoly.spad.pamphlet
+++ b/src/algebra/xlpoly.spad.pamphlet
@@ -1269,7 +1269,6 @@ r + r1 + r2
 --E 28
 )spool
  
-)spool
 )lisp (bye)
 @
 <<LiePolynomial.help>>=
diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html
index ca14fd6..34d60cf 100644
--- a/src/axiom-website/patches.html
+++ b/src/axiom-website/patches.html
@@ -791,6 +791,10 @@ regression file fixed created<br/>
 CATS hyperbolicrules.input added<br/>
 <a href="patches/20081209.01.tpd.patch">20081209.01.tpd.patch</a>
 bookvol10.3 add domains<br/>
+<a href="patches/20081210.01.tpd.patch">20081210.01.tpd.patch</a>
+bookvol10.3 add domains<br/>
+<a href="patches/20081211.01.tpd.patch">20081211.01.tpd.patch</a>
+<br/>
 
  </body>
 </html>
\ No newline at end of file
diff --git a/src/input/Makefile.pamphlet b/src/input/Makefile.pamphlet
index 45343ce..b341753 100644
--- a/src/input/Makefile.pamphlet
+++ b/src/input/Makefile.pamphlet
@@ -266,7 +266,7 @@ OUTS= ffrac.output     \
       tutchap2.output  tutchap3.output  tutchap4.output  \
       up.output        \
       vector.output    viewdef.output   \
-      wutset.output    xpbwpoly.output  \
+      wutset.output    \
       xpoly.output     xpr.output       \
       zdsolve.output   zimmer.output    zlindep.output 
 
@@ -687,7 +687,7 @@ FILES= ${OUT}/algaggr.input  ${OUT}/algbrbf.input    
${OUT}/algfacob.input \
        ${OUT}/vector.input   ${OUT}/vectors.input    ${OUT}/viewdef.input \
        ${OUT}/void.input     ${OUT}/wiggle.input   \
        ${OUT}/wutset.input \
-       ${OUT}/xpbwpoly.input ${OUT}/xpoly.input      ${OUT}/xpr.input \
+       ${OUT}/xpoly.input      ${OUT}/xpr.input \
        ${OUT}/zdsolve.input  ${OUT}/zimmer.input     ${OUT}/zlindep.input
 
 FILES2=${OUT}/arith.input    ${OUT}/bugs.input \
@@ -1043,7 +1043,7 @@ DOCFILES= \
   ${DOC}/vector.input.dvi      ${DOC}/vectors.input.dvi    \
   ${DOC}/viewdef.input.dvi     ${DOC}/void.input.dvi       \
   ${DOC}/wester.input.dvi      ${DOC}/wiggle.input.dvi     \
-  ${DOC}/wutset.input.dvi      ${DOC}/xpbwpoly.input.dvi   \
+  ${DOC}/wutset.input.dvi      \
   ${DOC}/xpoly.input.dvi       ${DOC}/xpr.input.dvi        \
   ${DOC}/zdsolve.input.dvi     ${DOC}/zimmer.input.dvi     \
   ${DOC}/zlindep.input.dvi  
diff --git a/src/input/cwmmt.input.pamphlet b/src/input/cwmmt.input.pamphlet
index c037652..2179513 100644
--- a/src/input/cwmmt.input.pamphlet
+++ b/src/input/cwmmt.input.pamphlet
@@ -98,6 +98,7 @@ CompWithMappingModeTest() : Exports == Implementation where
 \end{chunk}
 <<*>>=
 )spool cwmmt.output
+)sys cp $AXIOM/../../src/input/cwmmt.input.pamphlet .
 )lisp (tangle "cwmmt.input.pamphlet" "cwmmt.spad" "cwmmt.spad" )
 )co cwmmt.spad
 )set message test on
diff --git a/src/input/function.input.pamphlet 
b/src/input/function.input.pamphlet
index 8e2c181..52c3385 100644
--- a/src/input/function.input.pamphlet
+++ b/src/input/function.input.pamphlet
@@ -359,6 +359,8 @@ sinCosExpand(sin(x+y-2*z) * cos y)
 --R                                                     Type: Expression 
Integer
 --E 33
 
+)spool
+)lisp (bye)
 @
 \eject
 \begin{thebibliography}{99}
diff --git a/src/input/series.input.pamphlet b/src/input/series.input.pamphlet
index bd53d23..eb840eb 100644
--- a/src/input/series.input.pamphlet
+++ b/src/input/series.input.pamphlet
@@ -18,7 +18,8 @@
 )set message test on
 )set message auto off
 )clear all
---S 1 of 17
+
+@
 \section{Expression To Power Series}
 We compute series expansions of various functions using EXPR2UPS.
 
diff --git a/src/input/xpbwpoly.input.pamphlet 
b/src/input/xpbwpoly.input.pamphlet
deleted file mode 100644
index f838d24..0000000
--- a/src/input/xpbwpoly.input.pamphlet
+++ /dev/null
@@ -1,59 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/input XPBWPOLY.input}
-\author{The Axiom Team}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-<<*>>=
-)cl all
-
-a:Symbol := 'a
-b:Symbol := 'b
-RN    := Fraction(Integer)
-word   := OrderedFreeMonoid Symbol
-lword := LyndonWord(Symbol)
-base  := PoincareBirkhoffWittLyndonBasis Symbol
-dpoly := XDistributedPolynomial(Symbol, RN)
-rpoly := XRecursivePolynomial(Symbol, RN)
-lpoly := LiePolynomial(Symbol, RN)
-poly  := XPBWPolynomial(Symbol, RN)
-liste : List lword := LyndonWordsList([a,b], 6)
-0$poly
-1$poly
-p : poly := a
-q : poly := b
-pq: poly := p*q
-pq :: dpoly
-mirror pq
-ListOfTerms pq
-reductum pq
-leadingMonomial pq
-coefficients pq
-leadingTerm pq
-degree pq
-pq4:=exp(pq,4)
-log(pq4,4) - pq
-lp1 :lpoly := LiePoly liste.10
-lp2 :lpoly := LiePoly liste.11
-lp  :lpoly := [lp1, lp2]
-lpd1: dpoly := lp1
-lpd2: dpoly := lp2
-lpd : dpoly := lpd1 * lpd2 - lpd2 * lpd1
-lp :: dpoly - lpd
-p := 3 * lp
-q := lp1
-pq:= p * q
-pr:rpoly := p :: rpoly
-qr:rpoly := q :: rpoly
-pq :: rpoly - pr*qr
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}




reply via email to

[Prev in Thread] Current Thread [Next in Thread]