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