diff --git a/books/bookvol10.3.pamphlet b/books/bookvol10.3.pamphlet
index 2b35d0e..0015f59 100644
--- a/books/bookvol10.3.pamphlet
+++ b/books/bookvol10.3.pamphlet
@@ -119074,6 +119074,877 @@ SingletonAsOrderedSet(): OrderedSet with
\end{chunk}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{domain SEM SparseEchelonMatrix}
+
+\begin{chunk}{SparseEchelonMatrix.input}
+)set break resume
+)sys rm -f SparseEchelonMatrix.output
+)spool SparseEchelonMatrix.output
+)set message test on
+)set message auto off
+)clear all
+
+--S 1 of 1
+)show SparseEchelonMatrix
+--R SparseEchelonMatrix(C: OrderedSet,D: Ring) is a domain constructor
+--R Abbreviation for SparseEchelonMatrix is SEM
+--R This constructor is exposed in this frame.
+--R Issue )edit /research/ref/SEM.spad to see algebra source code for SEM
+--R
+--R------------------------------- Operations --------------------------------
+--R ?*? : (Matrix(D),%) -> % allIndices : % -> List(C)
+--R coerce : % -> Matrix(D) coerce : % -> OutputForm
+--R copy : % -> % deleteRow! : (%,Integer) -> Void
+--R elimZeroCols! : % -> Void elt : (%,Integer,C) -> D
+--R horizJoin : (%,%) -> % join : (%,%) -> %
+--R ncols : % -> NonNegativeInteger new : (List(C),Integer) -> %
+--R nrows : % -> NonNegativeInteger setelt! : (%,Integer,C,D) -> Void
+--R ?*? : (Matrix(Fraction(D)),%) -> % if D has INTDOM
+--R appendRow! : (%,Record(Indices: List(C),Entries: List(D))) -> Void
+--R consRow! : (%,Record(Indices: List(C),Entries: List(D))) -> Void
+--R extract : (%,Integer,Integer) -> %
+--R horizSplit : (%,C) -> Record(Left: %,Right: %)
+--R pivot : (%,Integer) -> Record(Index: C,Entry: D)
+--R pivots : % -> Record(Indices: List(C),Entries: List(D))
+--R primitiveRowEchelon : % -> Record(Ech: %,Lt: Matrix(Fraction(D)),Pivots: List(D),Rank: NonNegativeInteger) if D has GCDDOM
+--R purge! : (%,(C -> Boolean)) -> Void
+--R row : (%,Integer) -> Record(Indices: List(C),Entries: List(D))
+--R rowEchelon : % -> Record(Ech: %,Lt: Matrix(D),Pivots: List(D),Rank: NonNegativeInteger)
+--R setGcdMode : Symbol -> Symbol if D has GCDDOM
+--R setRow! : (%,Integer,List(C),List(D)) -> Void
+--R setRow! : (%,Integer,Record(Indices: List(C),Entries: List(D))) -> Void
+--R sortedPurge! : (%,(C -> Boolean)) -> Void
+--R
+--E 1
+
+)spool
+)lisp (bye)
+\end{chunk}
+\begin{chunk}{SparseEchelonMatrix.help}
+====================================================================
+SparseEchelonMatrix examples
+====================================================================
+
+See Also:
+o )show SparseEchelonMatrix
+
+\end{chunk}
+
+\pagehead{SparseEchelonMatrix}{SEM}
+\pagepic{ps/v103sparseechelonmatrix.eps}{SEM}{1.00}
+
+{\bf Exports:}\\
+\begin{tabular}{lll}
+\cross{SEM}{?*?} &
+\cross{SEM}{allIndices} &
+\cross{SEM}{coerce} \\
+\cross{SEM}{coerce} &
+\cross{SEM}{copy} &
+\cross{SEM}{deleteRow!} \\
+\cross{SEM}{elimZeroCols!} &
+\cross{SEM}{elt} &
+\cross{SEM}{horizJoin} \\
+\cross{SEM}{join} &
+\cross{SEM}{ncols} &
+\cross{SEM}{new} \\
+\cross{SEM}{nrows} &
+\cross{SEM}{setelt!} &
+\cross{SEM}{?*?} \\
+\cross{SEM}{appendRow!} &
+\cross{SEM}{consRow!} &
+\cross{SEM}{extract} \\
+\cross{SEM}{horizSplit} &
+\cross{SEM}{pivot} &
+\cross{SEM}{pivots} \\
+\cross{SEM}{primitiveRowEchelon} &
+\cross{SEM}{purge!} &
+\cross{SEM}{row} \\
+\cross{SEM}{rowEchelon} &
+\cross{SEM}{setGcdMode} &
+\cross{SEM}{setRow!} \\
+\cross{SEM}{setRow!} &
+\cross{SEM}{sortedPurge!}
+\end{tabular}
+
+\begin{chunk}{domain SEM SparseEchelonMatrix}
+)abbrev domain SEM SparseEchelonMatrix
+++ Description:
+++ \spad{SparseEchelonMatrix(C, D)} implements sparse matrices whose columns
+++ are enumerated by the \spadtype{OrderedSet} \spad{C} and whose entries
+++ belong to the \spadtype{GcdDomain} \spad{D}. The basic operation of
+++ this domain is the computation of an row echelon form. The used algorithm
+++ tries to maintain the sparsity and is especially adapted to matrices who
+++ are already close to a row echelon form.
+
+SparseEchelonMatrix(C : OrderedSet, D : Ring) : Cat == Def where
+
+ Sy ==> Symbol
+ L ==> List
+ V ==> Vector
+ VD ==> Vector D
+ MD ==> Matrix D
+ FD ==> Fraction D
+ MFD ==> Matrix FD
+ I ==> Integer
+ NNI ==> NonNegativeInteger
+ B ==> Boolean
+ OUT ==> OutputForm
+
+ ROWREC ==> Record(Indices : L C, Entries : L D)
+
+ iter ==> "iterated"::Sy
+ rand ==> "random"::Sy
+
+ Cat ==> CoercibleTo OUT with
+
+ shallowlyMutable
+ ++ Matrices may be altered destructively.
+
+ finiteAggregate
+ ++ Matrices are finite.
+
+ coerce : % -> MD
+ ++ \spad{coerce(A)} yields the matrix \spad{A} in the usual matrix type.
+
+ copy : % -> %
+ ++ \spad{copy(A)} returns a copy of the matrix \spad{A}.
+
+ ncols : % -> NNI
+ ++ \spad{ncols(A)} returns the number of columns of the matrix \spad{A}.
+
+ nrows : % -> NNI
+ ++ \spad{nrows(A)} returns the number of rows of the matrix \spad{A}.
+
+ allIndices : % -> L C
+ ++ \spad{allIndices(A)} returns all indices used for enumerating the
+ ++ columns of the matrix \spad{A}.
+
+ elimZeroCols! : % -> Void
+ ++ \spad{elimZeroCols!(A)} removes columns which contain only zeros.
+ ++ This effects basically only the value of \spad{allIndices(A)}.
+
+ purge! : (%, C-> B) -> Void
+ ++ \spad{purge!(A, crit)} eliminates all columns belonging to an index
+ ++ \spad{c} such that \spad{crit(c)} yields \spad{true}.
+
+ sortedPurge! : (%, C-> B) -> Void
+ ++ \spad{sortedPurge!(A, crit)} is like \spad{purge}, however, with the
+ ++ additional assumption that \spad{crit} respects the ordering of the
+ ++ indices.
+
+ new : (L C, I) -> %
+ ++ \spad{new(inds, nrows)} generates a new matrix with \spad{nrows}
+ ++ rows and columns enumerated by the indices \spad{inds}. The matrix
+ ++ is empty, i.e. the zero matrix.
+
+ elt : (%, I, C) -> D
+ ++ \spad{elt(A, i, c)} returns the entry of the matrix \spad{A} in row
+ ++ \spad{i} and in the column with index \spad{c}.
+
+ setelt! : (%, I, C, D) -> Void
+ ++ \spad{setelt!(A, i, c, d)} sets the entry of the matrix \spad{A} in
+ ++ row \spad{i} and in the column with index \spad{c} to the value
+ ++ \spad{d}.
+
+ row : (%, I) -> ROWREC
+ ++ \spad{row(A, i)} returns the \spad{i}-th row of the matrix \spad{A}.
+
+ setRow! : (%, I, ROWREC) -> Void
+ ++ \spad{setRow!(A, i, ind, ent)} sets the \spad{i}-th row of the matrix
+ ++ \spad{A} to the value \spad{r}.
+
+ setRow! : (%, I, L C, L D) -> Void
+ ++ \spad{setRow!(A, i, ind, ent)} sets the \spad{i}-th row of the matrix
+ ++ \spad{A}. Its indices are \spad{ind}; the entries \spad{ent}.
+
+ deleteRow! : (%, I) -> Void
+ ++ \spad{deleteRow(A, i)} deletes the \spad{i}-th row of the matrix
+ ++ \spad{A}.
+
+ consRow! : (%, ROWREC) -> Void
+ ++ \spad{consRow!(A, r)} inserts the row \spad{r} at the top of the
+ ++ matrix \spad{A}.
+
+ appendRow! : (%, ROWREC) -> Void
+ ++ \spad{appendRow!(A, r)} appends the row \spad{r} at the end of the
+ ++ matrix \spad{A}.
+
+ extract : (%, I, I) -> %
+ ++ \spad{extract(A, i1, i2)} extracts the rows \spad{i1} to \spad{i2}
+ ++ and returns them as a new matrix.
+
+ rowEchelon : % -> Record(Ech : %, Lt : MD, Pivots : L D, Rank : NNI)
+ ++ \spad{primitiveRowEchelon(A)} computes a row echelon form for the
+ ++ matrix \spad{A}. The algorithm used is fraction-free elimination.
+ ++ It is especially adapted to matrices already close to row echelon
+ ++ form. The transformation matrix, the used pivots and the rank of the
+ ++ matrix are also returned.
+
+ if D has GcdDomain then
+
+ setGcdMode : Sy -> Sy
+ ++ \spad{setGcdMode(s)} sets a new value for the flag deciding on
+ ++ the method used to compute gcd`s for lists. Possible values for
+ ++ \spad{s} are \spad{iterated} and \spad{random}.
+
+ primitiveRowEchelon: % -> Record(Ech:%, Lt:MFD, Pivots:L D, Rank:NNI)
+ ++ \spad{primitiveRowEchelon(A)} computes a row echelon form for the
+ ++ matrix \spad{A}. The algorithm used is fraction-free elimination.
+ ++ Every row is made primitive by division by the gcd. The algorithm
+ ++ is especially adapted to matrices already close to row echelon
+ ++ form. The transformation matrix, the used pivots and the rank of the
+ ++ matrix are also returned.
+
+ pivot : (%, I) -> Record(Index : C, Entry : D)
+ ++ \spad{pivot(A, i)} returns the leading entry of the \spad{i}-th row
+ ++ of the matrix \spad{A} together with its index.
+
+ pivots : % -> ROWREC
+ ++ \spad{pivots(A)} returns all leading entries of the matrix \spad{A}
+ ++ together with their indices.
+
+ "*" : (MD, %) -> %
+ ++ \spad{L*A} implements left multiplication with a usual matrix.
+
+ if D has IntegralDomain then
+
+ "*" : (MFD, %) -> %
+ ++ \spad{L*A} implements left multiplication with a usual matrix over
+ ++ the quotient field of \spad{D}.
+
+ join : (%, %) -> %
+ ++ \spad{join(A, B)} vertically concats the matrices \spad{A} and
+ ++ \spad{B}.
+
+ horizJoin : (%, %) -> %
+ ++ \spad{horizJoin(A, B)} horizontally concats the matrices \spad{A} and
+ ++ \spad{B}. It is assumed that all indices of \spad{B} are smaller than
+ ++ those of \spad{A}.
+
+ horizSplit : (%, C) -> Record(Left : %, Right : %)
+ ++ \spad{horizSplit(A, c)} splits the matrix \spad{A} into two at the
+ ++ column given by \spad{c}. The first column of the right matrix is
+ ++ enumerated by the first index less or equal to \spad{c}.
+
+ Def ==> add
+
+ minInd : I := minIndex([i for i in 1..1])
+ offset : I := minInd-1
+ emptyRec : ROWREC := [empty, empty]
+ noChecks? : B := D has lazyRepresentation -- flag for lazy representation
+ seed : I := 113 -- seed for random number generation
+ GCDmode : Sy := iter -- flag for gcd algorithm
+
+ greater(r1 : ROWREC, r2 : ROWREC) : B ==
+ empty? r1.Indices => false
+ empty? r2.Indices => true
+ first(r2.Indices) < first(r1.Indices)
+
+ -- -------------- --
+ -- Representation --
+ -- -------------- --
+
+ -- For efficiency reasons most checks for correct index ranges are omitted.
+
+ Rep := Record(NCols : NNI, NRows : NNI, AllInds : L C, Rows : V ROWREC)
+
+ ncols(A : %) : NNI == A.NCols
+
+ nrows(A : %) : NNI == A.NRows
+
+ allIndices(A : %) : L C == copy A.AllInds
+
+ row(A : %, i : I) : ROWREC ==
+ -- i < 0 or i > A.NRows => error "index out of range"
+ qelt(A.Rows, i)
+
+ setRow!(A : %, i : I, r : ROWREC) : Void ==
+ -- i < 0 or i > A.NRows => error "index out of range"
+ qsetelt!(A.Rows, i, r)
+ void
+
+ setRow!(A : %, i : I, inds : L C, ents : L D) : Void ==
+ -- i < 0 or i > A.NRows => error "index out of range"
+ -- #inds ^= #ents => error "improper row"
+ qsetelt!(A.Rows, i, [inds, ents])
+ void
+
+ new(inds : L C, n : I) : % ==
+ [#inds, n::NNI, inds, [copy emptyRec for i in 1..n]]
+
+ elt(A : %, i : I, c : C) : D ==
+ r := row(A, i)
+ pos := position(c, r.Indices)
+ pos < minInd => 0$D
+ qelt(r.Entries, pos)
+
+ setelt!(A : %, i : I, c : C, d : D) : Void ==
+ r := row(A, i)
+ pos := position(c, r.Indices)
+ if pos >= minInd then
+ qsetelt!(r.Entries, pos, d)
+ else
+ j := minInd
+ for ind in r.Indices while c < ind repeat
+ j := j+1
+ r.Indices := insert!(c, r.Indices, j)
+ r.Entries := insert!(d, r.Entries, j)
+ qsetelt!(A.Rows, i, r)
+ void
+
+ coerce(A : %) : MD ==
+ zero? A.NCols => error "cannot coerce matrix with zero columns"
+ AA : MD := new(A.NRows, A.NCols, 0$D)
+ for r in entries(A.Rows) for i in minRowIndex(AA).. repeat
+ inds := r.Indices
+ ents := r.Entries
+ for ind in A.AllInds for j in minColIndex(AA).. _
+ while not empty? inds repeat
+ if ind = first inds then
+ qsetelt!(AA, i, j, first ents)
+ inds := rest inds
+ ents := rest ents
+ AA
+
+ coerce(A : %) : OUT ==
+ zero? A.NCols => 0$D ::OUT
+ A::MD::OUT
+
+ copy(A : %) : % ==
+ resRows : V ROWREC := new(A.NRows, emptyRec)
+ for l in 1..A.NRows repeat
+ r := qelt(A.Rows, l)
+ qsetelt!(resRows, l, [copy r.Indices, copy r.Entries])
+ [A.NCols, A.NRows, copy A.AllInds, resRows]
+
+ -- ----------------------- --
+ -- Basic Matrix Operations --
+ -- ----------------------- --
+
+ elimZeroCols!(A : %) : Void ==
+ newInds : L C := empty
+ for r in entries(A.Rows) repeat
+ newInds := removeDuplicates! merge!((x, y) +-> y < x,
+ newInds, r.Indices)
+ A.AllInds := newInds
+ void
+
+ purge!(A : %, crit : C-> B) : Void ==
+ newInds : L C := empty
+ for c in A.AllInds repeat
+ if not crit c then
+ newInds := cons(c, newInds)
+ newInds := reverse! newInds
+ if #newInds ^= #A.AllInds then
+ A.AllInds := newInds
+ for l in 1..A.NRows repeat
+ r := qelt(A.Rows, l)
+ newInds : L C := empty
+ newEnts : L D := empty
+ for c in r.Indices for e in r.Entries repeat
+ if not crit c then
+ newInds := cons(c, newInds)
+ newEnts := cons(e, newEnts)
+ qsetelt!(A.Rows, l, [reverse! newInds, reverse! newEnts])
+ void
+
+ sortedPurge!(A : %, crit : C-> B) : Void ==
+ if crit first A.AllInds then
+ while not(empty? A.AllInds) and crit first A.AllInds repeat
+ A.AllInds := rest A.AllInds
+ for l in 1..A.NRows repeat
+ r := qelt(A.Rows, l)
+ while not(empty? r.Indices) and crit first r.Indices repeat
+ r.Indices := rest r.Indices
+ r.Entries := rest r.Entries
+ qsetelt!(A.Rows, l, r)
+ void
+
+ deleteRow!(A : %, i : I) : Void ==
+ i > A.NRows => A
+ nr := (A.NRows-1)::NNI
+ resRows : V ROWREC := new(nr, emptyRec)
+ for l in 1..(i-1) repeat
+ qsetelt!(resRows, l, qelt(A.Rows, l))
+ for l in (i+1)..A.NRows repeat
+ qsetelt!(resRows, l-1, qelt(A.Rows, l))
+ A.NRows := nr
+ A.Rows := resRows
+ void
+
+ consRow!(A : %, r : ROWREC) : Void ==
+ A.NRows := A.NRows + 1
+ newRows : L ROWREC := cons(r, entries A.Rows)
+ A.Rows := construct newRows
+ newInds := setDifference(r.Indices, A.AllInds)
+ if not empty? newInds then
+ A.AllInds := merge((x, y) +-> y < x, A.AllInds,
+ sort!((x, y) +-> y < x, newInds))
+ void
+
+ appendRow!(A : %, r : ROWREC) : Void ==
+ A.NRows := A.NRows + 1
+ newRows : L ROWREC := concat(entries A.Rows, r)
+ A.Rows := construct newRows
+ newInds := setDifference(r.Indices, A.AllInds)
+ if not empty? newInds then
+ A.AllInds := merge((x, y) +-> y < x, A.AllInds,
+ sort!((x, y) +-> y < x, newInds))
+ void
+
+ extract(A : %, i1 : I, i2 : I) : % ==
+ nr := (i2-i1+1)::NNI
+ resRows : V ROWREC := new(nr, emptyRec)
+ newInds : L C := empty
+ for i in i1..i2 repeat
+ qsetelt!(resRows, i-i1+1, row(A, i))
+ newInds := removeDuplicates! merge((x, y) +-> y < x,
+ newInds, row(A, i).Indices)
+ [A.NCols, nr, newInds, resRows]
+
+ join(A1 : %, A2 : %) : % ==
+ newInds := removeDuplicates! merge((x : C, y : C) : Boolean +-> y < x,
+ A1.AllInds, A2.AllInds)
+ newNRows := A1.NRows + A2.NRows
+ newRows : V ROWREC := new(newNRows, emptyRec)
+ for l in 1..A1.NRows repeat
+ qsetelt!(newRows, l, qelt(A1.Rows, l))
+ for l in 1..A2.NRows repeat
+ qsetelt!(newRows, A1.NRows+l, qelt(A2.Rows, l))
+ [#newInds, newNRows, newInds, newRows]
+
+ horizJoin(A1 : %, A2 : %) : % ==
+ A1.NRows ^= A2.NRows => error "incompatible dimensions in horizJoin"
+ newInds := append(A1.AllInds, A2.AllInds)
+ res : % := new(newInds, A1.NRows)
+ for i in 1..A1.NRows repeat
+ r1 := row(A1, i)
+ r2 := row(A2, i)
+ setRow!(res, i, append(r1.Indices, r2.Indices), _
+ append(r1.Entries, r2.Entries))
+ res
+
+ horizSplit(A : %, c : C) : Record(Left : %, Right : %) ==
+ rinds : L C := allIndices A
+ linds : L C := empty
+ while not(empty? rinds) and (first(rinds) > c) repeat
+ linds := cons(first(rinds), linds)
+ rinds := rest rinds
+ empty? linds => [new(linds, A.NRows), A]
+ linds := reverse! linds
+ empty? rinds => [A, new(rinds, A.NRows)]
+ LA : % := new(linds, A.NRows)
+ RA : % := new(rinds, A.NRows)
+ for i in 1..A.NRows repeat
+ r := row(A, i)
+ ri : L C := r.Indices
+ re : L D := r.Entries
+ li : L C := empty
+ le : L D := empty
+ while not(empty? ri) and (first(ri) > c) repeat
+ li := cons(first(ri), li)
+ le := cons(first re, le)
+ ri := rest ri
+ re := rest re
+ if not empty? li then
+ li := reverse! li
+ le := reverse! le
+ setRow!(LA, i, li, le)
+ if not empty? ri then
+ setRow!(RA, i, ri, re)
+ [LA, RA]
+
+ -- ----------- --
+ -- Row Echelon --
+ -- ----------- --
+
+ addRows(d1 : D, r1 : ROWREC, d2 : D, r2 : ROWREC) : ROWREC ==
+ -- Computes linear combination of two rows.
+ -- Local function.
+ empty? r1.Indices =>
+ one? d2 => r2
+ [r2.Indices, [d2*e2 for e2 in r2.Entries]]
+ empty? r2.Indices =>
+ one? d1 => r1
+ [r1.Indices, [d1*e1 for e1 in r1.Entries]]
+ resI : L C := empty
+ resE : L D := empty
+ lent1 : L D
+ lent2 : L D
+ if not(noChecks?) and one? d1 then
+ lent1 := r1.Entries
+ else
+ lent1 := [d1*e1 for e1 in r1.Entries]
+ if not(noChecks?) and one? d2 then
+ lent2 := copy r2.Entries
+ else
+ lent2 := [d2*e2 for e2 in r2.Entries]
+ lind2 := copy r2.Indices
+
+ for c1 in r1.Indices for e1 in lent1 repeat
+ while not(empty? lind2) and c1 < first(lind2) repeat
+ resI := cons(first lind2, resI)
+ resE := cons(first(lent2), resE)
+ lind2 := rest lind2
+ lent2 := rest lent2
+ if not(empty? lind2) and first(lind2) = c1 then
+ sum := e1+first(lent2)
+ if noChecks? or not zero? sum then
+ resI := cons(c1, resI)
+ resE := cons(sum, resE)
+ lind2 := rest lind2
+ lent2 := rest lent2
+ else
+ resI := cons(c1, resI)
+ resE := cons(e1, resE)
+
+ resI := concat!(reverse! resI, lind2)
+ resE := concat!(reverse! resE, lent2)
+ while not(empty? resE) and zero? first resE repeat
+ resI := rest resI
+ resE := rest resE
+ [resI, resE]
+
+ pivot(A : %, i : I) : Record(Index : C, Entry : D) ==
+ r := row(A, i)
+ empty? r.Indices => error "empty row"
+ [first r.Indices, first r.Entries]
+
+ pivots(A : %) : ROWREC ==
+ resI : L C := empty
+ resE : L D := empty
+ for r in entries A.Rows | not empty? r.Indices repeat
+ resI := cons(first r.Indices, resI)
+ resE := cons(first r.Entries, resE)
+ [reverse! resI, reverse! resE]
+
+ rowEchelon(AA : %) : Record(Ech : %, Lt : MD, Pivots : L D, Rank : NNI) ==
+ A := copy AA
+ LTr : MD := diagonalMatrix [1$D for i in 1..A.NRows]
+ Pivs : L D := empty
+
+ -- check pivots
+ for i in 1..A.NRows repeat
+ r := qelt(A.Rows, i)
+ changed? : B := false
+ while not(empty? r.Entries) and zero? first r.Entries repeat
+ r.Entries := rest r.Entries
+ r.Indices := rest r.Indices
+ changed? := true
+ if changed? then
+ qsetelt!(A.Rows, i, r)
+
+ -- sort rows by pivots (bubble sort)
+ sorted? : B := false
+ until sorted? repeat
+ sorted? := true
+ oldr := qelt(A.Rows, 1)
+ for i in 2..A.NRows repeat
+ newr := qelt(A.Rows, i)
+ if greater(newr, oldr) then
+ qsetelt!(A.Rows, i, oldr)
+ qsetelt!(A.Rows, i-1, newr)
+ swapRows!(LTr, i-1, i)
+ sorted? := false
+ else
+ oldr := newr
+
+ -- fraction-free elimination
+ finished? : B := false
+ pivlen, pivrow, rk : NNI
+ for i in 1..A.NRows until finished? repeat
+ r := qelt(A.Rows, i)
+ finished? := empty? r.Indices
+ if finished? then
+ rk : NNI := (i-1)::NNI
+ else -- search good pivot
+ pivind := first r.Indices
+ pivlen := #r.Indices
+ pivrow := i
+ k : I := 0
+ for j in (i+1)..A.NRows _
+ while not(empty? qelt(A.Rows, j).Indices) and _
+ pivind = first(qelt(A.Rows, j).Indices) repeat
+ len := #qelt(A.Rows, j).Indices
+ k := k+1
+ if len < pivlen then
+ pivlen := len
+ pivrow := j
+ piv : D := first qelt(A.Rows, pivrow).Entries
+ Pivs := cons(piv, Pivs)
+
+ -- elimination necessary?
+ if k > 0 then
+ if pivrow ^= i then
+ pr := qelt(A.Rows, pivrow)
+ qsetelt!(A.Rows, pivrow, qelt(A.Rows, i))
+ qsetelt!(A.Rows, i, pr)
+ swapRows!(LTr, i, pivrow)
+
+ -- elimination (and resorting of rows)
+ pr := copy qelt(A.Rows, i)
+ pr.Indices := rest pr.Indices
+ pr.Entries := rest pr.Entries
+ for j in (i+1)..(i+k) repeat
+ r := copy qelt(A.Rows, i+1)
+ c := first r.Entries
+ r.Indices := rest r.Indices
+ r.Entries := rest r.Entries
+ r := addRows(piv, r, -c, pr)
+ for l in 1..A.NRows repeat
+ f := piv*qelt(LTr, i+1, l) - c*qelt(LTr, i, l)
+ qsetelt!(LTr, i+1, l, f)
+ for l in (i+2)..(2*i+k+1-j) repeat
+ qsetelt!(A.Rows, l-1, qelt(A.Rows, l))
+ swapRows!(LTr, l-1, l)
+ for l in (2*i+k+2-j)..A.NRows _
+ while greater(qelt(A.Rows, l), r) repeat
+ qsetelt!(A.Rows, l-1, qelt(A.Rows, l))
+ swapRows!(LTr, l-1, l)
+ qsetelt!(A.Rows, l-1, r)
+
+ if not finished? then
+ rk : NNI := A.NRows
+ [A, LTr, Pivs, rk]
+
+ if D has GcdDomain then
+
+ setGcdMode(s : Sy) : Sy ==
+ tmp := GCDmode
+ (s = iter) or (s = rand) =>
+ GCDmode := s
+ tmp
+ error "unknown gcd mode"
+
+ randomGCD(le : L D) : D ==
+ -- Probabilistic technique.
+ #le = 2 => gcd(first le, second le)
+ f := first le
+ g := second le
+ l := rest rest le
+ while not empty? l repeat
+ one? first l => return 1$D
+ f := f + (1+random(113)$I)*first(l)
+ l := rest l
+ if not empty? l then
+ one? first l => return 1$D
+ g := g + (1+random(113)$I)*first(l)
+ l := rest l
+ h := gcd(f, g)
+ l := [h]
+ for e in le repeat
+ tmp := e exquo h
+ if tmp case "failed" then
+ l := cons(e, l)
+ one?(#l) => h
+ randomGCD l
+
+ iteratedGCD(le : L D) : D ==
+ -- Computes gcd iteratively
+ res := gcd(first le, second le)
+ l := rest rest le
+ while not(empty?(l) or one?(res)) repeat
+ res := gcd(res, first l)
+ l := rest l
+ res
+
+ makePrimitive(r : ROWREC) : Record(GCD : D, Row : ROWREC) ==
+ -- remove common gcd of row
+ le := r.Entries
+ one?(#le) => [first le, [r.Indices, [1$D]]]
+ g : D
+ if GCDmode = 'iterated then
+ g := iteratedGCD le
+ else
+ g := randomGCD le
+ one? g => [1, r]
+ le := [(e exquo g)::D for e in le]
+ [g, [r.Indices, le]]
+
+ primitiveRowEchelon(AA : %) : _
+ Record(Ech : %, Lt : MFD, Pivots : L D, Rank : NNI) ==
+ A := copy AA
+ LTr : MFD := diagonalMatrix [1$FD for i in 1..A.NRows]
+ Pivs : L D := empty
+
+ -- check pivots
+ for i in 1..A.NRows repeat
+ r := qelt(A.Rows, i)
+ changed? : B := false
+ while not(empty? r.Entries) and zero? first r.Entries repeat
+ r.Entries := rest r.Entries
+ r.Indices := rest r.Indices
+ changed? := true
+ if changed? then
+ qsetelt!(A.Rows, i, r)
+
+ -- sort rows by pivots (bubble sort)
+ sorted? : B := false
+ until sorted? repeat
+ sorted? := true
+ oldr := qelt(A.Rows, 1)
+ for i in 2..A.NRows repeat
+ newr := qelt(A.Rows, i)
+ if greater(newr, oldr) then
+ qsetelt!(A.Rows, i, oldr)
+ qsetelt!(A.Rows, i-1, newr)
+ swapRows!(LTr, i-1, i)
+ sorted? := false
+ else
+ oldr := newr
+
+ -- primitive fraction-free elimination
+ finished? : B := false
+ pivlen, pivrow, rk : NNI
+ for i in 1..A.NRows until finished? repeat
+ r := qelt(A.Rows, i)
+ finished? := empty? r.Indices
+ if finished? then
+ rk : NNI := (i-1)::NNI
+ else -- search good pivot
+ pivind := first r.Indices
+ pivlen := #r.Indices
+ pivrow := i
+ k : I := 0
+ for j in (i+1)..A.NRows _
+ while not(empty? qelt(A.Rows, j).Indices) and _
+ pivind = first(qelt(A.Rows, j).Indices) repeat
+ len := #qelt(A.Rows, j).Indices
+ k := k+1
+ if len < pivlen then
+ pivlen := len
+ pivrow := j
+
+ -- make row primitive
+ tmp := makePrimitive qelt(A.Rows, pivrow)
+ if not one? tmp.GCD then
+ qsetelt!(A.Rows, pivrow, tmp.Row)
+ q : FD := 1/tmp.GCD
+ for l in 1..A.NRows | not zero? qelt(LTr, pivrow, l) _
+ repeat
+ qsetelt!(LTr, pivrow, l, q*qelt(LTr, pivrow, l))
+ piv : D := first qelt(A.Rows, pivrow).Entries
+ Pivs := cons(piv, Pivs)
+
+ -- elimination necessary?
+ if k > 0 then
+ if pivrow ^= i then
+ pr := qelt(A.Rows, pivrow)
+ qsetelt!(A.Rows, pivrow, qelt(A.Rows, i))
+ qsetelt!(A.Rows, i, pr)
+ swapRows!(LTr, i, pivrow)
+
+ -- elimination (and resorting of rows)
+ pr := copy tmp.Row
+ pr.Indices := rest pr.Indices
+ pr.Entries := rest pr.Entries
+ for j in (i+1)..(i+k) repeat
+ r := copy qelt(A.Rows, i+1)
+ c := first r.Entries
+ r.Indices := rest r.Indices
+ r.Entries := rest r.Entries
+ r := addRows(piv, r, -c, pr)
+ for l in 1..A.NRows repeat
+ fd : FD := piv *$FD qelt(LTr, i+1, l) - _
+ (c*qelt(LTr, i, l))::FD
+ qsetelt!(LTr, i+1, l, fd)
+ for l in (i+2)..(2*i+k+1-j) repeat
+ qsetelt!(A.Rows, l-1, qelt(A.Rows, l))
+ swapRows!(LTr, l-1, l)
+ for l in (2*i+k+2-j)..A.NRows _
+ while greater(qelt(A.Rows, l), r) repeat
+ qsetelt!(A.Rows, l-1, qelt(A.Rows, l))
+ swapRows!(LTr, l-1, l)
+ qsetelt!(A.Rows, l-1, r)
+
+ if not finished? then
+ rk : NNI := A.NRows
+ [A, LTr, Pivs, rk]
+
+ -- -------------- --
+ -- Multiplication --
+ -- -------------- --
+
+ L : MD * AA : % ==
+ ncols(L) ^= AA.NRows => error "improper matrix dimensions"
+ A := copy AA
+ rlen := nrows L
+ res : % := new(A.AllInds, rlen)
+
+ for c in A.AllInds repeat
+ tmp : V D := new(rlen, 0$D)
+ for i in 1..A.NRows repeat
+ r := qelt(A.Rows, i)
+ inds := r.Indices
+ if not(empty? inds) and first(inds) = c then
+ for k in 1..rlen | not zero? qelt(L, k, i) repeat
+ qsetelt!(tmp, k, qelt(tmp, k) + qelt(L, k, i)* _
+ first(r.Entries))
+ r.Entries := rest r.Entries
+ r.Indices := rest inds
+ qsetelt!(A.Rows, i, r)
+ for k in 1..rlen | not zero? qelt(tmp, k) repeat
+ r := qelt(res.Rows, k)
+ r.Indices := cons(c, r.Indices)
+ r.Entries := cons(qelt(tmp, k), r.Entries)
+ qsetelt!(res.Rows, k, r)
+
+ for k in 1..rlen repeat
+ r := qelt(res.Rows, k)
+ r.Indices := reverse! r.Indices
+ r.Entries := reverse! r.Entries
+ qsetelt!(res.Rows, k, r)
+ res
+
+ if D has IntegralDomain then
+
+ mult(f : FD, d : D) : D ==
+ res := numer(f)*d
+ tmp := res exquo denom(f)
+ tmp case "failed" => error "cannot divide in mult"
+ tmp::D
+
+ L : MFD * AA : % ==
+ ncols(L) ^= AA.NRows => error "improper matrix dimensions"
+ A := copy AA
+ rlen := nrows L
+ res : % := new(A.AllInds, rlen)
+
+ for c in A.AllInds repeat
+ tmp : V FD := new(rlen, 0$FD)
+ for i in 1..A.NRows repeat
+ r := qelt(A.Rows, i)
+ inds := r.Indices
+ if not(empty? inds) and first(inds) = c then
+ for k in 1..rlen | not zero? qelt(L, k, i) repeat
+ qsetelt!(tmp, k, qelt(tmp, k) + qelt(L, k, i)* _
+ first(r.Entries))
+ r.Entries := rest r.Entries
+ r.Indices := rest inds
+ qsetelt!(A.Rows, i, r)
+ for k in 1..rlen | not zero? qelt(tmp, k) repeat
+ d : Union(D, "failed") := retractIfCan qelt(tmp, k)
+ d case "failed" => error "cannot divide in *"
+ r := qelt(res.Rows, k)
+ r.Indices := cons(c, r.Indices)
+ r.Entries := cons(d::D, r.Entries)
+ qsetelt!(res.Rows, k, r)
+
+ for k in 1..rlen repeat
+ r := qelt(res.Rows, k)
+ r.Indices := reverse! r.Indices
+ r.Entries := reverse! r.Entries
+ qsetelt!(res.Rows, k, r)
+ res
+
+\end{chunk}
+\begin{chunk}{SEM.dotabb}
+"SEM" [color="#88FF44",href="bookvol10.3.pdf#nameddest=SEM"]
+"ALIST" [color="#88FF44",href="bookvol10.3.pdf#nameddest=ALIST"]
+"SEM" -> "ALIST"
+
+\end{chunk}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{domain SMP SparseMultivariatePolynomial}
\begin{chunk}{SparseMultivariatePolynomial.input}
diff --git a/books/bookvol10.pamphlet b/books/bookvol10.pamphlet
index 53b0bea..c8f1684 100644
--- a/books/bookvol10.pamphlet
+++ b/books/bookvol10.pamphlet
@@ -11034,7 +11034,7 @@ LAYER16=\
${OUT}/RESULT.o ${OUT}/RFFACT.o ${OUT}/RMATRIX.o ${OUT}/ROMAN.o \
${OUT}/ROUTINE.o ${OUT}/RPOLCAT.o ${OUT}/RPOLCAT-.o ${OUT}/RULECOLD.o \
${OUT}/SAOS.o ${OUT}/SCELL.o \
- ${OUT}/SEGBIND.o ${OUT}/SET.o ${OUT}/SPECOUT.o \
+ ${OUT}/SEGBIND.o ${OUT}/SEM.o ${OUT}/SET.o ${OUT}/SPECOUT.o \
${OUT}/SQMATRIX.o ${OUT}/SWITCH.o ${OUT}/SYMS.o ${OUT}/SYMTAB.o \
${OUT}/SYSSOLP.o ${OUT}/UTSCAT.o ${OUT}/UTSCAT-.o ${OUT}/VARIABLE.o \
layer16done
@@ -12882,6 +12882,29 @@ LAYER16=\
/*"SEGBIND" -> {"A1AGG-"; "ISTRING"; "SRAGG-"; "FLAGG-"; "LNAGG-"}*/
/*"SEGBIND" -> "BOOLEAN"*/
+"SEM" [color="#88FF44",href="bookvol10.3.pdf#nameddest=SEM"]
+"SEM" -> "ALIST"
+/*"SEM" -> {"A1AGG", "A1AGG-", "ABELGRP", "ABELMON", "ABELMON-"; */
+/*"SEM" -> {ABELSG", "ABELSG-", "AGG", "AGG-", "ALGEBRA", "ALIST"; */
+/*"SEM" -> {BASTYPE", "BASTYPE-", "BMODULE", "BOOLEAN", "CABMON"; */
+/*"SEM" -> {CFCAT", "CHAR", "CHARNZ", "CHARZ", "CLAGG", "CLAGG-"; */
+/*"SEM" -> {COMRING", "DIFEXT", "DIFRING", "DIVRING", "ELAGG"; */
+/*"SEM" -> {ELAGG-", "ELTAB", "ELTAGG", "ELTAGG-", "ENTIRER"; */
+/*"SEM" -> {EUCDOM", "EVALAB", "FEVALAB", "FIELD", "FLAGG", "FLAGG-"; */
+/*"SEM" -> {FLINEXP", "FPATMAB", "GCDDOM", "HOAGG", "HOAGG-"; */
+/*"SEM" -> {IARRAY1", "IEVALAB", "ILIST", "INS", "INT", "INTDOM"; */
+/*"SEM" -> {ISTRING", "IVECTOR", "IXAGG", "IXAGG-", "KOERCE"; */
+/*"SEM" -> {KONVERT", "LINEXP", "LIST", "LMODULE", "LNAGG", "LNAGG-"; */
+/*"SEM" -> {LORER", "LSAGG", "LSAGG-", "MODULE", "MONOID", "MONOID-"; */
+/*"SEM" -> {NNI", "OAGROUP", "OAMON", "OASGP", "OCAMON", "OINTDOM"; */
+/*"SEM" -> {OM", "ORDRING", "ORDSET", "ORDSET-", "OUTFORM", "PATAB"; */
+/*"SEM" -> {PATMAB", "PDRING", "PFECAT", "PI", "PID", "PRIMARR"; */
+/*"SEM" -> {QFCAT", "RCAGG", "RCAGG-", "REAL", "REF", "RETRACT"; */
+/*"SEM" -> {RING", "RMODULE", "RNG", "SETCAT", "SETCAT-", "SGROUP"; */
+/*"SEM" -> {SGROUP-", "SINT", "SRAGG-", "STAGG", "STAGG-", "STEP"; */
+/*"SEM" -> {STRING", "SYMBOL", "TYPE", "UFD", "URAGG", "URAGG-"; */
+/*"SEM" -> {VECTCAT", "VECTCAT-", "VECTOR"; */
+
"SET" [color="#88FF44",href="bookvol10.3.pdf#nameddest=SET"]
/*"SET" -> {"FSAGG"; "DIAGG"; "DIOPS"; "BGAGG"; "HOAGG"; "AGG"; "TYPE"}*/
/*"SET" -> {"SETCAT"; "BASTYPE"; "KOERCE"; "EVALAB"; "IEVALAB"; "CLAGG"}*/
@@ -18935,6 +18958,7 @@ REGRESS= \
SmithNormalForm.regress \
SortedCache.regress \
SortPackage.regress \
+ SparseEchelonMatrix.regress \
SparseTable.regress \
SparseMultivariatePolynomial.regress \
SparseMultivariateTaylorSeries.regress \
diff --git a/books/bookvol5.pamphlet b/books/bookvol5.pamphlet
index ad1491b..b2d34a7 100644
--- a/books/bookvol5.pamphlet
+++ b/books/bookvol5.pamphlet
@@ -24569,6 +24569,7 @@ otherwise the new algebra won't be loaded by the interpreter when needed.
(|SimplifyAlgebraicNumberConvertPackage| . SIMPAN)
(|SingleInteger| . SINT)
(|SmithNormalForm| . SMITH)
+ (|SparseEchelonMatrix| . SEM)
(|SparseUnivariatePolynomialExpressions| . SUPEXPR)
(|SparseUnivariatePolynomialFunctions2| . SUP2)
(|SpecialOutputPackage| . SPECOUT)
diff --git a/books/ps/v103sparseechelonmatrix.eps b/books/ps/v103sparseechelonmatrix.eps
new file mode 100644
index 0000000..c1a38bf
--- /dev/null
+++ b/books/ps/v103sparseechelonmatrix.eps
@@ -0,0 +1,278 @@
+%!PS-Adobe-3.0 EPSF-3.0
+%%Creator: graphviz version 2.26.3 (20100126.1600)
+%%Title: pic
+%%Pages: 1
+%%BoundingBox: 36 36 100 152
+%%EndComments
+save
+%%BeginProlog
+/DotDict 200 dict def
+DotDict begin
+
+/setupLatin1 {
+mark
+/EncodingVector 256 array def
+ EncodingVector 0
+
+ISOLatin1Encoding 0 255 getinterval putinterval
+EncodingVector 45 /hyphen put
+
+% Set up ISO Latin 1 character encoding
+/starnetISO {
+ dup dup findfont dup length dict begin
+ { 1 index /FID ne { def }{ pop pop } ifelse
+ } forall
+ /Encoding EncodingVector def
+ currentdict end definefont
+} def
+/Times-Roman starnetISO def
+/Times-Italic starnetISO def
+/Times-Bold starnetISO def
+/Times-BoldItalic starnetISO def
+/Helvetica starnetISO def
+/Helvetica-Oblique starnetISO def
+/Helvetica-Bold starnetISO def
+/Helvetica-BoldOblique starnetISO def
+/Courier starnetISO def
+/Courier-Oblique starnetISO def
+/Courier-Bold starnetISO def
+/Courier-BoldOblique starnetISO def
+cleartomark
+} bind def
+
+%%BeginResource: procset graphviz 0 0
+/coord-font-family /Times-Roman def
+/default-font-family /Times-Roman def
+/coordfont coord-font-family findfont 8 scalefont def
+
+/InvScaleFactor 1.0 def
+/set_scale {
+ dup 1 exch div /InvScaleFactor exch def
+ scale
+} bind def
+
+% styles
+/solid { [] 0 setdash } bind def
+/dashed { [9 InvScaleFactor mul dup ] 0 setdash } bind def
+/dotted { [1 InvScaleFactor mul 6 InvScaleFactor mul] 0 setdash } bind def
+/invis {/fill {newpath} def /stroke {newpath} def /show {pop newpath} def} bind def
+/bold { 2 setlinewidth } bind def
+/filled { } bind def
+/unfilled { } bind def
+/rounded { } bind def
+/diagonals { } bind def
+
+% hooks for setting color
+/nodecolor { sethsbcolor } bind def
+/edgecolor { sethsbcolor } bind def
+/graphcolor { sethsbcolor } bind def
+/nopcolor {pop pop pop} bind def
+
+/beginpage { % i j npages
+ /npages exch def
+ /j exch def
+ /i exch def
+ /str 10 string def
+ npages 1 gt {
+ gsave
+ coordfont setfont
+ 0 0 moveto
+ (\() show i str cvs show (,) show j str cvs show (\)) show
+ grestore
+ } if
+} bind def
+
+/set_font {
+ findfont exch
+ scalefont setfont
+} def
+
+% draw text fitted to its expected width
+/alignedtext { % width text
+ /text exch def
+ /width exch def
+ gsave
+ width 0 gt {
+ [] 0 setdash
+ text stringwidth pop width exch sub text length div 0 text ashow
+ } if
+ grestore
+} def
+
+/boxprim { % xcorner ycorner xsize ysize
+ 4 2 roll
+ moveto
+ 2 copy
+ exch 0 rlineto
+ 0 exch rlineto
+ pop neg 0 rlineto
+ closepath
+} bind def
+
+/ellipse_path {
+ /ry exch def
+ /rx exch def
+ /y exch def
+ /x exch def
+ matrix currentmatrix
+ newpath
+ x y translate
+ rx ry scale
+ 0 0 1 0 360 arc
+ setmatrix
+} bind def
+
+/endpage { showpage } bind def
+/showpage { } def
+
+/layercolorseq
+ [ % layer color sequence - darkest to lightest
+ [0 0 0]
+ [.2 .8 .8]
+ [.4 .8 .8]
+ [.6 .8 .8]
+ [.8 .8 .8]
+ ]
+def
+
+/layerlen layercolorseq length def
+
+/setlayer {/maxlayer exch def /curlayer exch def
+ layercolorseq curlayer 1 sub layerlen mod get
+ aload pop sethsbcolor
+ /nodecolor {nopcolor} def
+ /edgecolor {nopcolor} def
+ /graphcolor {nopcolor} def
+} bind def
+
+/onlayer { curlayer ne {invis} if } def
+
+/onlayers {
+ /myupper exch def
+ /mylower exch def
+ curlayer mylower lt
+ curlayer myupper gt
+ or
+ {invis} if
+} def
+
+/curlayer 0 def
+
+%%EndResource
+%%EndProlog
+%%BeginSetup
+14 default-font-family set_font
+1 setmiterlimit
+% /arrowlength 10 def
+% /arrowwidth 5 def
+
+% make sure pdfmark is harmless for PS-interpreters other than Distiller
+/pdfmark where {pop} {userdict /pdfmark /cleartomark load put} ifelse
+% make '<<' and '>>' safe on PS Level 1 devices
+/languagelevel where {pop languagelevel}{1} ifelse
+2 lt {
+ userdict (<<) cvn ([) cvn load put
+ userdict (>>) cvn ([) cvn load put
+} if
+
+%%EndSetup
+setupLatin1
+%%Page: 1 1
+%%PageBoundingBox: 36 36 100 152
+%%PageOrientation: Portrait
+0 0 1 beginpage
+gsave
+36 36 64 116 boxprim clip newpath
+1 1 set_scale 0 rotate 40 41 translate
+0.16355 0.45339 0.92549 graphcolor
+newpath -4 -5 moveto
+-4 112 lineto
+61 112 lineto
+61 -5 lineto
+closepath fill
+1 setlinewidth
+0.16355 0.45339 0.92549 graphcolor
+newpath -4 -5 moveto
+-4 112 lineto
+61 112 lineto
+61 -5 lineto
+closepath stroke
+% SEM
+gsave
+[ /Rect [ 1 72 55 108 ]
+ /Border [ 0 0 0 ]
+ /Action << /Subtype /URI /URI (bookvol10.3.pdf#nameddest=SEM) >>
+ /Subtype /Link
+/ANN pdfmark
+0.27273 0.73333 1 nodecolor
+newpath 55 108 moveto
+1 108 lineto
+1 72 lineto
+55 72 lineto
+closepath fill
+1 setlinewidth
+filled
+0.27273 0.73333 1 nodecolor
+newpath 55 108 moveto
+1 108 lineto
+1 72 lineto
+55 72 lineto
+closepath stroke
+0 0 0 nodecolor
+14 /Times-Roman set_font
+12.5 86.4 moveto 31 (SEM) alignedtext
+grestore
+% ALIST
+gsave
+[ /Rect [ 0 0 56 36 ]
+ /Border [ 0 0 0 ]
+ /Action << /Subtype /URI /URI (bookvol10.3.pdf#nameddest=ALIST) >>
+ /Subtype /Link
+/ANN pdfmark
+0.27273 0.73333 1 nodecolor
+newpath 56 36 moveto
+0 36 lineto
+0 0 lineto
+56 0 lineto
+closepath fill
+1 setlinewidth
+filled
+0.27273 0.73333 1 nodecolor
+newpath 56 36 moveto
+0 36 lineto
+0 0 lineto
+56 0 lineto
+closepath stroke
+0 0 0 nodecolor
+14 /Times-Roman set_font
+7.5 14.4 moveto 41 (ALIST) alignedtext
+grestore
+% SEM->ALIST
+gsave
+1 setlinewidth
+0 0 0 edgecolor
+newpath 28 71.83 moveto
+28 64.13 28 54.97 28 46.42 curveto
+stroke
+0 0 0 edgecolor
+newpath 31.5 46.41 moveto
+28 36.41 lineto
+24.5 46.41 lineto
+closepath fill
+1 setlinewidth
+solid
+0 0 0 edgecolor
+newpath 31.5 46.41 moveto
+28 36.41 lineto
+24.5 46.41 lineto
+closepath stroke
+grestore
+endpage
+showpage
+grestore
+%%PageTrailer
+%%EndPage: 1
+%%Trailer
+end
+restore
+%%EOF
diff --git a/changelog b/changelog
index 6bf0843..a90755b 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,14 @@
+20140906 tpd src/axiom-website/patches.html 20140906.01.tpd.patch
+20140906 tpd books/bookvol10.3 add SparseEchelonMatrix domain
+20140906 tpd books/bookvol10 add SparseEchelonMatrix domain
+20140906 tpd books/bookvol5 add SparseEchelonMatrix domain
+20140906 tpd books/ps/v103sparseechelonmatrix.eps added
+20140906 tpd src/share/algebra/browse.daase add SparseEchelonMatrix
+20140906 tpd src/share/algebra/category.daase add SparseEchelonMatrix
+20140906 tpd src/share/algebra/dependents.daase/index.kaf add SEM
+20140906 tpd src/share/algebra/interp.daase add SparseEchelonMatrix
+20140906 tpd src/share/algebra/operation.daase add SparseEchelonMatrix
+20140906 tpd src/share/algebra/users.daase/index.kaf add SparseEchelonMatrix
20140904 tpd src/axiom-website/patches.html 20140904.01.tpd.patch
20140904 tpd books/ps/v104graphviz.eps
20140904 tpd books/bookvol10.3
diff --git a/patch b/patch
index dbdf69b..f45c72c 100644
--- a/patch
+++ b/patch
@@ -1,3 +1,3 @@
-books/bookvol10.4 add graphviz package
+books/bookvol10.3 add SparseEchelonMatrix domain
-Add a package to create encapsulated postscript (.eps) graphs
+Add the SparseEchelonMatrix (SEM)
diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html
index f195772..3ea8770 100644
--- a/src/axiom-website/patches.html
+++ b/src/axiom-website/patches.html
@@ -4628,6 +4628,8 @@ src/input/Makefile typo fix
src/axiom-website/videos.html add Building Test Cases video
20140904.01.tpd.patch
books/bookvol10.4 add graphviz package
+20140906.01.tpd.patch
+books/bookvol10.3 add SparseEchelonMatrix domain