axiom-developer
[Top][All Lists]

## [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}
+\pagepic{ps/v103euclideanmodularring.ps}{EMR}{1.00}
+\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.
+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
+
+
+    --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
+          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
+             yv:=reduction(invlcy*yv,xm)
+           dy:=degree yv
+           xv:=x.val
+           while (d:=degree xv - dy)>=0 repeat
+                 xv:=reduction(fmecg(xv,d::NonNegativeInteger,
+                 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
+             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
+        (leadingCoefficient(x.val) = 1) => x
+        invlcx * x
+
+      unitNormal x ==
+--        zero?(x) or one?(leadingCoefficient(x.val)) => [1, x, 1]
+        zero?(x) or ((leadingCoefficient(x.val)) = 1) => [1, x, 1]
+        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}
\pagepic{ps/v103exponentialexpansion.ps}{EXPEXPAN}{1.00}
@@ -26690,6 +26801,82 @@ GeneralDistributedMultivariatePolynomial(vl,R,E):
public == private where

@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain GMODPOL GeneralModulePolynomial}
+\pagepic{ps/v103generalmodulepolynomial.ps}{GMODPOL}{1.00}
+\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
+        leadingExponent: $-> E + ++ leadingExponent(x) \undocumented + leadingIndex:$ -> IS
+        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 ==>
+
+      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}
+\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.
+    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
@
\pagepic{ps/v103linearordinarydifferentialoperator.ps}{LODO}{1.00}
+\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
@
\pagepic{ps/v103linearordinarydifferentialoperator2.ps}{LODO2}{1.00}
+\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>>=
+)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
+
+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
+    --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}
+\pagepic{ps/v103modularfield.ps}{MODFIELD}{1.00}
+\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.
+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}
+\pagepic{ps/v103modularring.ps}{MODRING}{1.00}
+\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.
+
+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
+
+    --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
+        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}
+\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}
+    scale: (%,F) -> %
+      ++ scale(m,h) returns \spad{scale(h) * m}
+    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)}
+    eval: (%,P1F) -> P1F
+      ++ eval(m,x) returns \spad{(a*x + b)/(c*x + d)}
+      ++ where \spad{m = moebius(a,b,c,d)}
+
+
+    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}
+\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:
+++  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
+
+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(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.
+
+        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>>=
+)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
+
+  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
+
+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}
+        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
+        remove: ( S -> Boolean ,%,Integer) -> %
+          ++ remove(p,ms,number) removes at most \spad{number} copies of
+          ++ if \spad{number} is positive, all of them if
+          ++ \spad{number} equals zero, and all but at most \spad{-number} if
+        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
+        remove_!: ( S -> Boolean ,%,Integer) -> %
+          ++ remove!(p,ms,number) removes destructively at most \spad{number}
+          ++ copies of elements x such that \spad{p(x)} is
+          ++ \spad{number} equals zero, and all but at most \spad{-number} if
+
+
+        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
+
+@
+\pagepic{ps/v103multivariatepolynomial.ps}{MPOLY}{1.00}
+\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>>=
+)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
+--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
+--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.
+
+        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.
+
+    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
+
+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
+      ++ 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
+          univariate(red p,mvar)
+
+--  a local functions, used in next definition
+      unlikeUnivReconstruct(u:SparseUnivariatePolynomial(%),mvar:VarSet):% ==
+        zero? (d:=degree u) => coefficient(u,0)
+            unlikeUnivReconstruct(reductum u,mvar)
+
+      multivariate(u:SparseUnivariatePolynomial(%),mvar:VarSet):% ==
+        ground? u => coefficient(u,0)
+        uu:=u
+        while not zero? uu repeat
+          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"
+
+--         p case R => p
+
+      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
+           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}
+\pagepic{ps/v103sparsemultivariatetaylorseries.ps}{SMTS}{1.00}
+\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) -> %
+        ++ to v and having c as the constant of integration.
+        ++ The evaluation of \spad{f()} is delayed.
+
+
+    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)
+
+      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
+        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
+          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
+
+@
+\pagepic{ps/v103squarematrix.ps}{SQMATRIX}{1.00}
+\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 ->$
+      ++ 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
+--    ++ 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}
+\pagepic{ps/v103taylorseries.ps}{TS}{1.00}
+\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) -> %
+        ++ to v and having c as the constant of integration.
+        ++ The evaluation of \spad{f()} is delayed.
+
+    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>>=
@@ -45407,6 +50643,7 @@ Note that this code is not included in the generated
<<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
<<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
<<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
<<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
<<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
<<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
<<domain QUEUE Queue>>

+<<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
<<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>>
index 4dbec57..f18549a 100644
@@ -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
-coefficients 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}