Initial population.

1 files changed, 2731 insertions, 0 deletions

diff --git a/src/booklets/Sorting.booklet b/src/booklets/Sorting.booklet
new file mode 100755
index 00000000..2f43240f
--- /dev/null
+++ b/src/booklets/Sorting.booklet
@@ -0,0 +1,2731 @@
+\documentclass{article}
+\usepackage{/home/axiomgnu/new/mnt/linux/bin/tex/noweb}
+\begin{document}
+\title{Sorting Facilities}
+\author{Timothy Daly}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+This is a survey of the explicitly mentioned sorting algorithms
+in the Axiom algebra code. Note that there are cases of "embedded"
+sorts as in the {\bf chainSubResultants} method from the
+{\bf PseudoRemainderSequence \cite{1}} domain. There are also
+implicit sorts as items are added to lists individually in sorted
+order.
+\subsection{aggcat.spad}
+\subsubsection{FiniteLinearAggregate}
+++ A finite linear aggregate is a linear aggregate of finite length.
+++ The finite property of the aggregate adds several exports to the
+++ list of exports from {\bf LinearAggregate} such as
+++ {\bf reverse}, {\bf sort}, and so on.
+<<FiniteLinearAggregate>>=
+FiniteLinearAggregate(S:Type): Category == LinearAggregate S with
+   finiteAggregate
+   sort: ((S,S)->B,%) -> %
+      ++ sort(p,a) returns a copy of {\bf a} sorted 
+      ++ using total ordering predicate p.
+   if S has OrderedSet then
+      OrderedSet
+      sort: % -> %
+	++ sort(u) returns an u with elements in ascending order.
+	++ Note: {\bf sort(u) = sort(<=,u)}.
+   if % has shallowlyMutable then
+      sort_!: ((S,S)->B,%) -> %
+	++ sort!(p,u) returns u with its elements ordered by p.
+      if S has OrderedSet then sort_!: % -> %
+	++ sort!(u) returns u with its elements in ascending order.
+ add
+    if S has OrderedSet then
+      sort l	  == sort(_<$S, l)
+
+    if % has shallowlyMutable then
+      sort(f, l) == sort_!(f, copy l)
+      if S has OrderedSet then
+	sort_! l == sort_!(_<$S, l)
+
+@
+\subsubsection{OneDimensionalArrayAggregate}
+One-dimensional-array aggregates serves as models for one-dimensional arrays.
+Categorically, these aggregates are finite linear aggregates
+with the {\bf shallowlyMutable} property, that is, any component of
+the array may be changed without affecting the
+identity of the overall array.
+Array data structures are typically represented by a fixed area in storage and
+therefore cannot efficiently grow or shrink on demand as can list structures
+(see however {\bf FlexibleArray} for a data structure which
+is a cross between a list and an array).
+Iteration over, and access to, elements of arrays is extremely fast
+(and often can be optimized to open-code).
+Insertion and deletion however is generally slow since an entirely new
+data structure must be created for the result.
+
+<<OneDimensionalArrayAggregate>>=
+OneDimensionalArrayAggregate(S:Type): Category ==
+    FiniteLinearAggregate S with shallowlyMutable
+  add
+    sort_!(f, a) == quickSort(f, a)$FiniteLinearAggregateSort(S, %)
+
+@
+\subsubsection{ListAggregate}
+A list aggregate is a model for a linked list data structure.
+A linked list is a versatile
+data structure. Insertion and deletion are efficient and
+searching is a linear operation.
+
+{\bf ListAggregate} uses the operation {\bf split\_!} which is inherited from
+{\bf StreamAggregate} which inherites it from it's implementing Domain
+{\bf UnaryRecursiveAggregate}. 
+A unary-recursive aggregate is a one where nodes may have either
+0 or 1 children.
+This aggregate models, though not precisely, a linked
+list possibly with a single cycle.
+A node with one children models a non-empty list, with the
+{\bf value} of the list designating the head, or {\bf first}, of the
+list, and the child designating the tail, or {\bf rest}, of the list.
+A node with no child then designates the empty list.
+Since these aggregates are recursive aggregates, they may be cyclic.
+
+There we see the definition of {\bf split\_!} as:
+
+<<UnaryRecursiveAggregate>>=
+UnaryRecursiveAggregate(S:Type): Category == RecursiveAggregate S with
+   rest: % -> %
+      ++ rest(u) returns an aggregate consisting of all but the first
+      ++ element of u
+      ++ (equivalently, the next node of u).
+   rest: (%,N) -> %
+      ++ rest(u,n) returns the \axiom{n}th (n >= 0) node of u.
+      ++ Note: \axiom{rest(u,0) = u}.
+   if % has shallowlyMutable then
+      split_!: (%,I) -> %
+	 ++ split!(u,n) splits u into two aggregates: \axiom{v = rest(u,n)}
+	 ++ and \axiom{w = first(u,n)}, returning \axiom{v}.
+	 ++ Note: afterwards \axiom{rest(u,n)} returns \axiom{empty()}.
+ add
+  rest(x, n) ==
+    for i in 1..n repeat
+      empty? x => error "Index out of range"
+      x := rest x
+    x
+
+  if S has SetCategory then
+    split_!(p, n) ==
+      n < 1 => error "index out of range"
+      p := rest(p, (n - 1)::N)
+      q := rest p
+      setrest_!(p, empty())
+      q
+
+
+@
+<<ListAggregate>>=
+ListAggregate(S:Type): Category == Join(StreamAggregate S,
+   FiniteLinearAggregate S, ExtensibleLinearAggregate S) with
+ add
+   mergeSort: ((S, S) -> B, %, I) -> %
+   mergeSort(f, p, n) ==
+     if n = 2 and f(first rest p, first p) then p := reverse_! p
+     n < 3 => p
+     l := (n quo 2)::N
+     q := split_!(p, l)
+     p := mergeSort(f, p, l)
+     q := mergeSort(f, q, n - l)
+     merge_!(f, p, q)
+
+   reverse_! x ==
+     empty? x => x
+     empty?(y := rest x) => x
+     setrest_!(x, empty())
+     while not empty? y repeat
+       z := rest y
+       setrest_!(y, x)
+       x := y
+       y := z
+     x
+
+   merge_!(f, p, q) ==
+     empty? p => q
+     empty? q => p
+     eq?(p, q) => error "cannot merge a list into itself"
+     if f(first p, first q)
+       then (r := t := p; p := rest p)
+       else (r := t := q; q := rest q)
+     while not empty? p and not empty? q repeat
+       if f(first p, first q)
+	 then (setrest_!(t, p); t := p; p := rest p)
+	 else (setrest_!(t, q); t := q; q := rest q)
+     setrest_!(t, if empty? p then q else p)
+     r
+
+   sort_!(f, l)	      == mergeSort(f, l, #l)
+
+   list x		   == concat(x, empty())
+   reduce(f, x)		   ==
+     empty? x => error "reducing over an empty list needs the 3 argument form"
+     reduce(f, rest x, first x)
+   merge(f, p, q)	   == merge_!(f, copy p, copy q)
+
+   select_!(f, x) ==
+     while not empty? x and not f first x repeat x := rest x
+     empty? x => x
+     y := x
+     z := rest y
+     while not empty? z repeat
+       if f first z then (y := z; z := rest z)
+		    else (z := rest z; setrest_!(y, z))
+     x
+
+   insert_!(s:S, x:%, i:I) ==
+     i < (m := minIndex x) => error "index out of range"
+     i = m => concat(s, x)
+     y := rest(x, (i - 1 - m)::N)
+     z := rest y
+     setrest_!(y, concat(s, z))
+     x
+
+   insert_!(w:%, x:%, i:I) ==
+     i < (m := minIndex x) => error "index out of range"
+     i = m => concat_!(w, x)
+     y := rest(x, (i - 1 - m)::N)
+     z := rest y
+     setrest_!(y, w)
+     concat_!(y, z)
+     x
+
+   remove_!(f:S -> B, x:%) ==
+     while not empty? x and f first x repeat x := rest x
+     empty? x => x
+     p := x
+     q := rest x
+     while not empty? q repeat
+       if f first q then q := setrest_!(p, rest q)
+		    else (p := q; q := rest q)
+     x
+
+   delete_!(x:%, i:I) ==
+     i < (m := minIndex x) => error "index out of range"
+     i = m => rest x
+     y := rest(x, (i - 1 - m)::N)
+     setrest_!(y, rest(y, 2))
+     x
+
+   delete_!(x:%, i:U) ==
+     (l := lo i) < (m := minIndex x) => error "index out of range"
+     h := if hasHi i then hi i else maxIndex x
+     h < l => x
+     l = m => rest(x, (h + 1 - m)::N)
+     t := rest(x, (l - 1 - m)::N)
+     setrest_!(t, rest(t, (h - l + 2)::N))
+     x
+
+   find(f, x) ==
+     while not empty? x and not f first x repeat x := rest x
+     empty? x => "failed"
+     first x
+
+   position(f:S -> B, x:%) ==
+     for k in minIndex(x).. while not empty? x and not f first x repeat
+       x := rest x
+     empty? x => minIndex(x) - 1
+     k
+
+   sorted?(f, l) ==
+     empty? l => true
+     p := rest l
+     while not empty? p repeat
+       not f(first l, first p) => return false
+       p := rest(l := p)
+     true
+
+   reduce(f, x, i) ==
+     r := i
+     while not empty? x repeat (r := f(r, first x); x := rest x)
+     r
+
+   if S has SetCategory then
+      reduce(f, x, i,a) ==
+	r := i
+	while not empty? x and r ^= a repeat
+	  r := f(r, first x)
+	  x := rest x
+	r
+
+   new(n, s) ==
+     l := empty()
+     for k in 1..n repeat l := concat(s, l)
+     l
+
+   map(f, x, y) ==
+     z := empty()
+     while not empty? x and not empty? y repeat
+       z := concat(f(first x, first y), z)
+       x := rest x
+       y := rest y
+     reverse_! z
+
+   copy x ==
+     y := empty()
+     for k in 0.. while not empty? x repeat
+       k = cycleMax and cyclic? x => error "cyclic list"
+       y := concat(first x, y)
+       x := rest x
+     reverse_! y
+
+   copyInto_!(y, x, s) ==
+     s < (m := minIndex y) => error "index out of range"
+     z := rest(y, (s - m)::N)
+     while not empty? z and not empty? x repeat
+       setfirst_!(z, first x)
+       x := rest x
+       z := rest z
+     y
+
+   if S has SetCategory then
+     position(w, x, s) ==
+       s < (m := minIndex x) => error "index out of range"
+       x := rest(x, (s - m)::N)
+       for k in s.. while not empty? x and w ^= first x repeat
+	 x := rest x
+       empty? x => minIndex x - 1
+       k
+
+     removeDuplicates_! l ==
+       p := l
+       while not empty? p repeat
+	 p := setrest_!(p, remove_!(#1 = first p, rest p))
+       l
+
+   if S has OrderedSet then
+     x < y ==
+	while not empty? x and not empty? y repeat
+	  first x ^= first y => return(first x < first y)
+	  x := rest x
+	  y := rest y
+	empty? x => not empty? y
+	false
+
+
+@
+\subsection{alql.spad}
+\subsubsection{DataList}
+This domain provides some nice functions on lists
+<<DataList>>=
+DataList(S:OrderedSet) : Exports == Implementation where
+  Exports == ListAggregate(S) with
+    elt: (%,"sort") -> %
+      ++ {\bf l.sort} returns {\bf l} with elements sorted.
+      ++ Note: {\bf l.sort = sort(l)}
+  Implementation == List(S) add
+    elt(x,"sort") == sort(x)
+
+@
+\subsection{carten.spad}
+\subsubsection{CartesianTensor}
+This is an instance of an "embedded sort" that should probably
+call one of the general purpose sort routines. 
+
+{\bf CartesianTensor(minix,dim,R)} provides Cartesian tensors with
+components belonging to a commutative ring R.  These tensors
+can have any number of indices.  Each index takes values from
+{\bf minix} to {\bf minix + dim - 1}.
+<<CartesianTensor>>=
+CartesianTensor(minix, dim, R): Exports == Implementation where
+    Exports ==> Join(GradedAlgebra(R, NNI), GradedModule(I, NNI)) with
+    Implementation ==> add
+        INDEX ==> Vector Integer  -- 1-based entries from minix..minix+dim-1
+
+        -- permsign!(v) = 1, 0, or -1  according as
+        -- v is an even, is not, or is an odd permutation of minix..minix+#v-1.
+        permsign_!(v: INDEX): Integer ==
+            -- sum minix..minix+#v-1.
+            maxix := minix+#v-1
+            psum  := (((maxix+1)*maxix - minix*(minix-1)) exquo 2)::Integer
+            -- +/v ^= psum => 0
+            n := 0
+            for i in 1..#v repeat n := n + v.i
+            n ^= psum => 0
+            -- Bubble sort!  This is pretty grotesque.
+            totTrans: Integer := 0
+            nTrans:   Integer := 1
+            while nTrans ^= 0 repeat
+                nTrans := 0
+                for i in 1..#v-1 for j in 2..#v repeat
+                    if v.i > v.j then
+                        nTrans := nTrans + 1
+                        e := v.i; v.i := v.j; v.j := e
+                totTrans := totTrans + nTrans
+            for i in 1..dim repeat
+                if v.i ^= minix+i-1 then return 0
+            odd? totTrans => -1
+            1
+
+@
+\subsection{clifford.spad}
+\subsubsection{CliffordAlgebra}
+Examples of {\bf Clifford Algebras} are: gaussians, quaternions, exterior 
+algebras and spin algebras.
+<<CliffordAlgebra>>=
+CliffordAlgebra(n, K, Q): T == Impl where
+    n: PositiveInteger
+    K: Field
+    Q: QuadraticForm(n, K)
+ 
+    PI ==> PositiveInteger
+    NNI==> NonNegativeInteger
+ 
+    T ==> Join(Ring, Algebra(K), VectorSpace(K)) with
+        e: PI -> %
+          ++ e(n) produces the appropriate unit element.
+        monomial: (K, List PI) -> %
+          ++ monomial(c,[i1,i2,...,iN]) produces the value given by
+          ++ \spad{c*e(i1)*e(i2)*...*e(iN)}.
+        coefficient:  (%, List PI) -> K
+          ++ coefficient(x,[i1,i2,...,iN])  extracts the coefficient of
+          ++ \spad{e(i1)*e(i2)*...*e(iN)} in x.
+        recip: % -> Union(%, "failed")
+          ++ recip(x) computes the multiplicative inverse of x or "failed"
+          ++ if x is not invertible.
+ 
+    Impl ==> add
+        Qeelist :=  [Q unitVector(i::PositiveInteger) for i in 1..n]
+        dim     :=  2**n
+ 
+        Rep     := PrimitiveArray K
+ 
+        New     ==> new(dim, 0$K)$Rep
+ 
+        x, y, z: %
+        c: K
+        m: Integer
+ 
+        characteristic() == characteristic()$K
+        dimension()      == dim::CardinalNumber
+ 
+        x = y ==
+            for i in 0..dim-1 repeat
+                if x.i ^= y.i then return false
+            true
+ 
+        x + y == (z := New; for i in 0..dim-1 repeat z.i := x.i + y.i; z)
+        x - y == (z := New; for i in 0..dim-1 repeat z.i := x.i - y.i; z)
+        - x   == (z := New; for i in 0..dim-1 repeat z.i := - x.i; z)
+        m * x == (z := New; for i in 0..dim-1 repeat z.i := m*x.i; z)
+        c * x == (z := New; for i in 0..dim-1 repeat z.i := c*x.i; z)
+ 
+        0            == New
+        1            == (z := New; z.0 := 1; z)
+        coerce(m): % == (z := New; z.0 := m::K; z)
+        coerce(c): % == (z := New; z.0 := c; z)
+ 
+        e b ==
+            b::NNI > n => error "No such basis element"
+            iz := 2**((b-1)::NNI)
+            z := New; z.iz := 1; z
+ 
+        -- The ei*ej products could instead be precomputed in
+        -- a (2**n)**2 multiplication table.
+        addMonomProd(c1: K, b1: NNI, c2: K, b2: NNI, z: %): % ==
+            c  := c1 * c2
+            bz := b2
+            for i in 0..n-1 | bit?(b1,i) repeat
+                -- Apply rule  ei*ej = -ej*ei for i^=j
+                k := 0
+                for j in i+1..n-1 | bit?(b1, j) repeat k := k+1
+                for j in 0..i-1   | bit?(bz, j) repeat k := k+1
+                if odd? k then c := -c
+                -- Apply rule  ei**2 = Q(ei)
+                if bit?(bz,i) then
+                    c := c * Qeelist.(i+1)
+                    bz:= (bz - 2**i)::NNI
+                else
+                    bz:= bz + 2**i
+            z.bz := z.bz + c
+            z
+ 
+        x * y ==
+            z := New
+            for ix in 0..dim-1 repeat
+                if x.ix ^= 0 then for iy in 0..dim-1 repeat
+                    if y.iy ^= 0 then addMonomProd(x.ix,ix,y.iy,iy,z)
+            z
+ 
+        canonMonom(c: K, lb: List PI): Record(coef: K, basel: NNI) ==
+            -- 0. Check input
+            for b in lb repeat b > n => error "No such basis element"
+ 
+            -- 1. Apply identity ei*ej = -ej*ei, i^=j.
+            -- The Rep assumes n is small so bubble sort is ok.
+            -- Using bubble sort keeps the exchange info obvious.
+            wasordered   := false
+            exchanges := 0
+            while not wasordered repeat
+                wasordered := true
+                for i in 1..#lb-1 repeat
+                    if lb.i > lb.(i+1) then
+                        t := lb.i; lb.i := lb.(i+1); lb.(i+1) := t
+                        exchanges := exchanges + 1
+                        wasordered := false
+            if odd? exchanges then c := -c
+ 
+            -- 2. Prepare the basis element
+            -- Apply identity ei*ei = Q(ei).
+            bz := 0
+            for b in lb repeat
+                bn := (b-1)::NNI
+                if bit?(bz, bn) then
+                    c := c * Qeelist bn
+                    bz:= ( bz - 2**bn )::NNI
+                else
+                    bz:= bz + 2**bn
+            [c, bz::NNI]
+ 
+        monomial(c, lb) ==
+            r := canonMonom(c, lb)
+            z := New
+            z r.basel := r.coef
+            z
+        coefficient(z, lb) ==
+            r := canonMonom(1, lb)
+            r.coef = 0 => error "Cannot take coef of 0"
+            z r.basel/r.coef
+ 
+        Ex ==> OutputForm
+ 
+        coerceMonom(c: K, b: NNI): Ex ==
+            b = 0 => c::Ex
+            ml := [sub("e"::Ex, i::Ex) for i in 1..n | bit?(b,i-1)]
+            be := reduce("*", ml)
+            c = 1 => be
+            c::Ex * be
+        coerce(x): Ex ==
+            tl := [coerceMonom(x.i,i) for i in 0..dim-1 | x.i^=0]
+            null tl => "0"::Ex
+            reduce("+", tl)
+
+
+	localPowerSets(j:NNI): List(List(PI)) ==
+	  l: List List PI := list []
+	  j = 0 => l
+	  Sm := localPowerSets((j-1)::NNI)
+          Sn: List List PI := []
+	  for x in Sm repeat Sn := cons(cons(j pretend PI, x),Sn)
+	  append(Sn, Sm)
+
+	powerSets(j:NNI):List List PI == map(reverse, localPowerSets j)
+
+	Pn:List List PI := powerSets(n)
+
+	recip(x: %): Union(%, "failed") ==
+	  one:% := 1
+	  -- tmp:c := x*yC - 1$C
+	  rhsEqs : List K := []
+	  lhsEqs: List List K := []
+	  lhsEqi: List K
+	  for pi in Pn repeat
+	    rhsEqs := cons(coefficient(one, pi), rhsEqs)
+
+	    lhsEqi := []
+	    for pj in Pn repeat
+		lhsEqi := cons(coefficient(x*monomial(1,pj),pi),lhsEqi)
+	    lhsEqs := cons(reverse(lhsEqi),lhsEqs)
+	  ans := particularSolution(matrix(lhsEqs),
+            vector(rhsEqs))$LinearSystemMatrixPackage(K, Vector K, Vector K, Matrix K)
+          ans case "failed" => "failed"
+	  ansP := parts(ans)
+	  ansC:% := 0
+	  for pj in Pn repeat
+	    cj:= first ansP
+	    ansP := rest ansP
+	    ansC := ansC + cj*monomial(1,pj)
+	  ansC
+
+@
+\subsection{defaults.spad}
+\subsubsection{FiniteLinearAggregateSort}
+++  This package exports 3 sorting algorithms which work over 
+++  FiniteLinearAggregates.
+<<FiniteLinearAggregateSort>>=
+FiniteLinearAggregateSort(S, V): Exports == Implementation where
+  S: Type
+  V: FiniteLinearAggregate(S) with shallowlyMutable
+ 
+  B ==> Boolean
+  I ==> Integer
+ 
+  Exports ==> with
+    quickSort: ((S, S) -> B, V) -> V
+      ++ quickSort(f, agg) sorts the aggregate agg with the ordering function
+      ++ f using the quicksort algorithm.
+    heapSort : ((S, S) -> B, V) -> V
+      ++ heapSort(f, agg) sorts the aggregate agg with the ordering function
+      ++ f using the heapsort algorithm.
+    shellSort: ((S, S) -> B, V) -> V
+      ++ shellSort(f, agg) sorts the aggregate agg with the ordering function
+      ++ f using the shellSort algorithm.
+ 
+  Implementation ==> add
+    siftUp   : ((S, S) -> B, V, I, I) -> Void
+    partition: ((S, S) -> B, V, I, I, I) -> I
+    QuickSort: ((S, S) -> B, V, I, I) -> V
+ 
+    quickSort(l, r) == QuickSort(l, r, minIndex r, maxIndex r)
+ 
+    siftUp(l, r, i, n) ==
+      t := qelt(r, i)
+      while (j := 2*i+1) < n repeat
+        if (k := j+1) < n and l(qelt(r, j), qelt(r, k)) then j := k
+        if l(t,qelt(r,j)) then
+           qsetelt_!(r, i, qelt(r, j))
+           qsetelt_!(r, j, t)
+           i := j
+        else leave
+ 
+    heapSort(l, r) ==
+      not zero? minIndex r => error "not implemented"
+      n := (#r)::I
+      for k in shift(n,-1) - 1 .. 0 by -1 repeat siftUp(l, r, k, n)
+      for k in n-1 .. 1 by -1 repeat
+         swap_!(r, 0, k)
+         siftUp(l, r, 0, k)
+      r
+ 
+    partition(l, r, i, j, k) ==
+      -- partition r[i..j] such that r.s <= r.k <= r.t
+      x := qelt(r, k)
+      t := qelt(r, i)
+      qsetelt_!(r, k, qelt(r, j))
+      while i < j repeat
+         if l(x,t) then
+           qsetelt_!(r, j, t)
+           j := j-1
+           t := qsetelt_!(r, i, qelt(r, j))
+         else (i := i+1; t := qelt(r, i))
+      qsetelt_!(r, j, x)
+      j
+ 
+    QuickSort(l, r, i, j) ==
+      n := j - i
+      if one? n and l(qelt(r, j), qelt(r, i)) then swap_!(r, i, j)
+      n < 2 => return r
+      -- for the moment split at the middle item
+      k := partition(l, r, i, j, i + shift(n,-1))
+      QuickSort(l, r, i, k - 1)
+      QuickSort(l, r, k + 1, j)
+ 
+    shellSort(l, r) ==
+      m := minIndex r
+      n := maxIndex r
+      -- use Knuths gap sequence: 1,4,13,40,121,...
+      g := 1
+      while g <= (n-m) repeat g := 3*g+1
+      g := g quo 3
+      while g > 0 repeat
+         for i in m+g..n repeat
+            j := i-g
+            while j >= m and l(qelt(r, j+g), qelt(r, j)) repeat
+               swap_!(r,j,j+g)
+               j := j-g
+         g := g quo 3
+      r
+ 
+
+@
+\subsection{e04agents.spad}
+\subsubsection{e04AgentsPackage}
+<<e04AgentsPackage>>=
+e04AgentsPackage(): E == I where
+
+++ Author: Brian Dupee
+++ Date Created: February 1996
+++ Date Last Updated: June 1996
+++ Basic Operations: simple? linear?, quadratic?, nonLinear?
+++ Description:
+++ \axiomType{e04AgentsPackage} is a package of numerical agents to be used
+++ to investigate attributes of an input function so as to decide the
+++ \axiomFun{measure} of an appropriate numerical optimization routine.
+
+  E ==> with
+
+    finiteBound:(LOCDF,DF) -> LDF 
+      ++ finiteBound(l,b) repaces all instances of an infinite entry in
+      ++ \axiom{l} by a finite entry \axiom{b} or \axiom{-b}.
+    sortConstraints:NOA -> NOA
+      ++ sortConstraints(args) uses a simple bubblesort on the list of
+      ++ constraints using the degree of the expression on which to sort.
+      ++ Of course, it must match the bounds to the constraints.
+    sumOfSquares:EDF -> Union(EDF,"failed")
+      ++ sumOfSquares(f) returns either an expression for which the square is
+      ++ the original function of "failed".
+    splitLinear:EDF -> EDF 
+      ++ splitLinear(f) splits the linear part from an expression which it
+      ++ returns.
+    simpleBounds?:LEDF -> Boolean
+      ++ simpleBounds?(l) returns true if the list of expressions l are
+      ++ simple.
+    linear?:LEDF -> Boolean
+      ++ linear?(l) returns true if all the bounds l are either linear or
+      ++ simple.
+    linear?:EDF -> Boolean
+      ++ linear?(e) tests if \axiom{e} is a linear function.
+    linearMatrix:(LEDF, NNI) -> MDF
+      ++ linearMatrix(l,n) returns a matrix of coefficients of the linear
+      ++ functions in \axiom{l}.  If l is empty, the matrix has at least one
+      ++ row.
+    linearPart:LEDF -> LEDF
+      ++ linearPart(l) returns the list of linear functions of \axiom{l}.
+    nonLinearPart:LEDF -> LEDF
+      ++ nonLinearPart(l) returns the list of non-linear functions of \axiom{l}.
+    quadratic?:EDF -> Boolean
+      ++ quadratic?(e) tests if \axiom{e} is a quadratic function.
+    variables:LSA -> LS
+      ++ variables(args) returns the list of variables in \axiom{args.lfn}
+    varList:(EDF,NNI) -> LS
+      ++ varList(e,n) returns a list of \axiom{n} indexed variables with name
+      ++ as in \axiom{e}.
+    changeNameToObjf:(Symbol,Result) -> Result
+      ++ changeNameToObjf(s,r) changes the name of item \axiom{s} in \axiom{r}
+      ++ to objf.
+    expenseOfEvaluation:LSA -> F
+      ++ expenseOfEvaluation(o) returns the intensity value of the 
+      ++ cost of evaluating the input set of functions.  This is in terms 
+      ++ of the number of ``operational units''.  It returns a value 
+      ++ in the range [0,1].
+    optAttributes:Union(noa:NOA,lsa:LSA) -> List String
+      ++ optAttributes(o) is a function for supplying a list of attributes
+      ++ of an optimization problem.
+
+  I ==> add
+
+    import ExpertSystemToolsPackage, ExpertSystemContinuityPackage
+
+    sumOfSquares2:EFI -> Union(EFI,"failed")
+    nonLinear?:EDF -> Boolean
+    finiteBound2:(OCDF,DF) -> DF 
+    functionType:EDF -> String
+
+    finiteBound2(a:OCDF,b:DF):DF ==
+      not finite?(a) =>
+        positive?(a) => b
+        -b
+      retract(a)@DF
+
+    finiteBound(l:LOCDF,b:DF):LDF == [finiteBound2(i,b) for i in l]
+
+    sortConstraints(args:NOA):NOA ==
+      Args := copy args
+      c:LEDF := Args.cf
+      l:LOCDF := Args.lb
+      u:LOCDF := Args.ub
+      m:INT := (# c) - 1      
+      n:INT := (# l) - m
+      for j in m..1 by -1 repeat
+        for i in 1..j repeat
+          s:EDF := c.i
+          t:EDF := c.(i+1)
+          if linear?(t) and (nonLinear?(s) or quadratic?(s)) then
+            swap!(c,i,i+1)$LEDF
+            swap!(l,n+i-1,n+i)$LOCDF
+            swap!(u,n+i-1,n+i)$LOCDF
+      Args
+        
+    changeNameToObjf(s:Symbol,r:Result):Result ==
+      a := remove!(s,r)$Result
+      a case Any =>
+        insert!([objf@Symbol,a],r)$Result
+        r
+      r
+
+    sum(a:EDF,b:EDF):EDF == a+b
+
+    variables(args:LSA): LS == variables(reduce(sum,(args.lfn)))
+
+    sumOfSquares(f:EDF):Union(EDF,"failed") ==
+      e := edf2efi(f)
+      s:Union(EFI,"failed") := sumOfSquares2(e)
+      s case EFI =>
+        map(fi2df,s)$EF2(FI,DF)
+      "failed"
+
+    sumOfSquares2(f:EFI):Union(EFI,"failed") ==
+      p := retractIfCan(f)@Union(PFI,"failed")
+      p case PFI => 
+        r := squareFreePart(p)$PFI
+        (p=r)@Boolean => "failed"
+        tp := totalDegree(p)$PFI
+        tr := totalDegree(r)$PFI
+        t := tp quo tr
+        found := false
+        q := r
+        for i in 2..t by 2 repeat
+          s := q**2
+          (s=p)@Boolean => 
+            found := true
+            leave
+          q := r**i
+        if found then 
+          q :: EFI
+        else
+          "failed"
+      "failed"
+
+    splitLinear(f:EDF):EDF ==
+      out := 0$EDF
+      (l := isPlus(f)$EDF) case LEDF =>
+        for i in l repeat
+          if not quadratic? i then
+            out := out + i
+        out
+      out
+
+    edf2pdf(f:EDF):PDF == (retract(f)@PDF)$EDF
+
+    varList(e:EDF,n:NNI):LS ==
+      s := name(first(variables(edf2pdf(e))$PDF)$LS)$Symbol
+      [subscript(s,[t::OutputForm]) for t in expand([1..n])$Segment(Integer)]
+
+    functionType(f:EDF):String ==
+      n := #(variables(f))$EDF
+      p := (retractIfCan(f)@Union(PDF,"failed"))$EDF
+      p case PDF =>
+        d := totalDegree(p)$PDF
+        one?(n*d) => "simple"
+        one?(d) => "linear"
+        (d=2)@Boolean => "quadratic"
+        "non-linear"
+      "non-linear"
+     
+    simpleBounds?(l: LEDF):Boolean ==
+      a := true
+      for e in l repeat
+        not (functionType(e) = "simple")@Boolean => 
+          a := false
+          leave
+      a
+
+    simple?(e:EDF):Boolean == (functionType(e) = "simple")@Boolean
+
+    linear?(e:EDF):Boolean == (functionType(e) = "linear")@Boolean
+
+    quadratic?(e:EDF):Boolean == (functionType(e) = "quadratic")@Boolean
+
+    nonLinear?(e:EDF):Boolean == (functionType(e) = "non-linear")@Boolean
+
+    linear?(l: LEDF):Boolean ==
+      a := true
+      for e in l repeat
+        s := functionType(e)
+        (s = "quadratic")@Boolean or (s = "non-linear")@Boolean => 
+          a := false
+          leave
+      a
+
+    simplePart(l:LEDF):LEDF == [i for i in l | simple?(i)]
+
+    linearPart(l:LEDF):LEDF == [i for i in l | linear?(i)]
+
+    nonLinearPart(l:LEDF):LEDF ==
+      [i for i in l | not linear?(i) and not simple?(i)]
+
+    linearMatrix(l:LEDF, n:NNI):MDF ==
+      empty?(l) => mat([],n)
+      L := linearPart l
+      M := zero(max(1,# L)$NNI,n)$MDF
+      vars := varList(first(l)$LEDF,n)
+      row:INT := 1
+      for a in L repeat
+        for j in monomials(edf2pdf(a))$PDF repeat
+          col:INT := 1
+          for c in vars repeat
+            if ((first(variables(j)$PDF)$LS)=c)@Boolean then
+              M(row,col):= first(coefficients(j)$PDF)$LDF
+            col := col+1
+        row := row + 1
+      M
+
+    expenseOfEvaluation(o:LSA):F ==
+      expenseOfEvaluation(vector(copy o.lfn)$VEDF)
+
+    optAttributes(o:Union(noa:NOA,lsa:LSA)):List String ==
+      o case noa =>
+        n := o.noa
+        s1:String := "The object function is " functionType(n.fn)
+        if empty?(n.lb) then
+          s2:String := "There are no bounds on the variables" 
+        else
+          s2:String := "There are simple bounds on the variables"
+        c := n.cf
+        if empty?(c) then
+          s3:String := "There are no constraint functions"
+        else
+          t := #(c)
+          lin := #(linearPart(c))
+          nonlin := #(nonLinearPart(c))
+          s3:String := "There are " string(lin)$String " linear and "_
+                          string(nonlin)$String " non-linear constraints"
+        [s1,s2,s3]
+      l := o.lsa
+      s:String := "non-linear"
+      if linear?(l.lfn) then
+        s := "linear"
+      ["The object functions are " s]
+
+@
+\subsection{expexpan.spad}
+++   UnivariatePuiseuxSeriesWithExponentialSingularity is a domain used to
+++   represent functions with essential singularities.  Objects in this
+++   domain are sums, where each term in the sum is a univariate Puiseux
+++   series times the exponential of a univariate Puiseux series.  Thus,
+++   the elements of this domain are sums of expressions of the form
+++   \spad{g(x) * exp(f(x))}, where g(x) is a univariate Puiseux series
+++   and f(x) is a univariate Puiseux series with no terms of non-negative
+++   degree.
+\subsubsection{UnivariatePuiseuxSeriesWithExponentialSingularity}
+<<UnivariatePuiseuxSeriesWithExponentialSingularity>>=
+UnivariatePuiseuxSeriesWithExponentialSingularity(R,FE,var,cen):_
+  Exports == Implementation where
+  R   : Join(OrderedSet,RetractableTo Integer,_
+             LinearlyExplicitRingOver Integer,GcdDomain)
+  FE  : Join(AlgebraicallyClosedField,TranscendentalFunctionCategory,_
+             FunctionSpace R)
+  var : Symbol
+  cen : FE
+  B       ==> Boolean
+  I       ==> Integer
+  L       ==> List
+  RN      ==> Fraction Integer
+  UPXS    ==> UnivariatePuiseuxSeries(FE,var,cen)
+  EXPUPXS ==> ExponentialOfUnivariatePuiseuxSeries(FE,var,cen)
+  OFE     ==> OrderedCompletion FE
+  Result  ==> Union(OFE,"failed")
+  PxRec   ==> Record(k: Fraction Integer,c:FE)
+  Term    ==> Record(%coef:UPXS,%expon:EXPUPXS,%expTerms:List PxRec)
+    -- the %expTerms field is used to record the list of the terms (a 'term'
+    -- records an exponent and a coefficient) in the exponent %expon
+  TypedTerm ==> Record(%term:Term,%type:String)
+    -- a term together with a String which tells whether it has an infinite,
+    -- zero, or unknown limit as var -> cen+
+  TRec    ==> Record(%zeroTerms: List Term,_
+                     %infiniteTerms: List Term,_
+                     %failedTerms: List Term,_
+                     %puiseuxSeries: UPXS)
+  SIGNEF  ==> ElementaryFunctionSign(R,FE)
+
+  Exports ==> Join(FiniteAbelianMonoidRing(UPXS,EXPUPXS),IntegralDomain) with
+    limitPlus : % -> Union(OFE,"failed")
+      ++ limitPlus(f(var)) returns \spad{limit(var -> cen+,f(var))}.
+    dominantTerm : % -> Union(TypedTerm,"failed")
+      ++ dominantTerm(f(var)) returns the term that dominates the limiting
+      ++ behavior of \spad{f(var)} as \spad{var -> cen+} together with a
+      ++ \spadtype{String} which briefly describes that behavior.  The
+      ++ value of the \spadtype{String} will be \spad{"zero"} (resp.
+      ++ \spad{"infinity"}) if the term tends to zero (resp. infinity)
+      ++ exponentially and will \spad{"series"} if the term is a
+      ++ Puiseux series.
+
+  Implementation ==> PolynomialRing(UPXS,EXPUPXS) add
+    makeTerm : (UPXS,EXPUPXS) -> Term
+    coeff : Term -> UPXS
+    exponent : Term -> EXPUPXS
+    exponentTerms : Term -> List PxRec
+    setExponentTerms_! : (Term,List PxRec) -> List PxRec
+    computeExponentTerms_! : Term -> List PxRec
+    terms : % -> List Term
+    sortAndDiscardTerms: List Term -> TRec
+    termsWithExtremeLeadingCoef : (L Term,RN,I) -> Union(L Term,"failed")
+    filterByOrder: (L Term,(RN,RN) -> B) -> Record(%list:L Term,%order:RN)
+    dominantTermOnList : (L Term,RN,I) -> Union(Term,"failed")
+    iDominantTerm : L Term -> Union(Record(%term:Term,%type:String),"failed")
+
+    retractIfCan f ==
+      (numberOfMonomials f = 1) and (zero? degree f) => leadingCoefficient f
+      "failed"
+
+    recip f ==
+      numberOfMonomials f = 1 =>
+        monomial(inv leadingCoefficient f,- degree f)
+      "failed"
+
+    makeTerm(coef,expon) == [coef,expon,empty()]
+    coeff term == term.%coef
+    exponent term == term.%expon
+    exponentTerms term == term.%expTerms
+    setExponentTerms_!(term,list) == term.%expTerms := list
+    computeExponentTerms_! term ==
+      setExponentTerms_!(term,entries complete terms exponent term)
+
+    terms f ==
+      -- terms with a higher order singularity will appear closer to the
+      -- beginning of the list because of the ordering in EXPPUPXS;
+      -- no "expnonent terms" are computed by this function
+      zero? f => empty()
+      concat(makeTerm(leadingCoefficient f,degree f),terms reductum f)
+
+    sortAndDiscardTerms termList ==
+      -- 'termList' is the list of terms of some function f(var), ordered
+      -- so that terms with a higher order singularity occur at the
+      -- beginning of the list.
+      -- This function returns lists of candidates for the "dominant
+      -- term" in 'termList', i.e. the term which describes the
+      -- asymptotic behavior of f(var) as var -> cen+.
+      -- 'zeroTerms' will contain terms which tend to zero exponentially
+      -- and contains only those terms with the lowest order singularity.
+      -- 'zeroTerms' will be non-empty only when there are no terms of
+      -- infinite or series type.
+      -- 'infiniteTerms' will contain terms which tend to infinity
+      -- exponentially and contains only those terms with the highest
+      -- order singularity.
+      -- 'failedTerms' will contain terms which have an exponential
+      -- singularity, where we cannot say whether the limiting value
+      -- is zero or infinity. Only terms with a higher order sigularity
+      -- than the terms on 'infiniteList' are included.
+      -- 'pSeries' will be a Puiseux series representing a term without an
+      -- exponential singularity.  'pSeries' will be non-zero only when no
+      -- other terms are known to tend to infinity exponentially
+      zeroTerms : List Term := empty()
+      infiniteTerms : List Term := empty()
+      failedTerms : List Term := empty()
+      -- we keep track of whether or not we've found an infinite term
+      -- if so, 'infTermOrd' will be set to a negative value
+      infTermOrd : RN := 0
+      -- we keep track of whether or not we've found a zero term
+      -- if so, 'zeroTermOrd' will be set to a negative value
+      zeroTermOrd : RN := 0
+      ord : RN := 0; pSeries : UPXS := 0  -- dummy values
+      while not empty? termList repeat
+        -- 'expon' is a Puiseux series
+        expon := exponent(term := first termList)
+        -- quit if there is an infinite term with a higher order singularity
+        (ord := order(expon,0)) > infTermOrd => leave "infinite term dominates"
+        -- if ord = 0, we've hit the end of the list
+        (ord = 0) =>
+          -- since we have a series term, don't bother with zero terms
+          leave(pSeries := coeff(term); zeroTerms := empty())
+        coef := coefficient(expon,ord)
+        -- if we can't tell if the lowest order coefficient is positive or
+        -- negative, we have a "failed term"
+        (signum := sign(coef)$SIGNEF) case "failed" =>
+          failedTerms := concat(term,failedTerms)
+          termList := rest termList
+        -- if the lowest order coefficient is positive, we have an
+        -- "infinite term"
+        (sig := signum :: Integer) = 1 =>
+          infTermOrd := ord
+          infiniteTerms := concat(term,infiniteTerms)
+          -- since we have an infinite term, don't bother with zero terms
+          zeroTerms := empty()
+          termList := rest termList
+        -- if the lowest order coefficient is negative, we have a
+        -- "zero term" if there are no infinite terms and no failed
+        -- terms, add the term to 'zeroTerms'
+        if empty? infiniteTerms then
+          zeroTerms :=
+            ord = zeroTermOrd => concat(term,zeroTerms)
+            zeroTermOrd := ord
+            list term
+        termList := rest termList
+      -- reverse "failed terms" so that higher order singularities
+      -- appear at the beginning of the list
+      [zeroTerms,infiniteTerms,reverse_! failedTerms,pSeries]
+
+    termsWithExtremeLeadingCoef(termList,ord,signum) ==
+      -- 'termList' consists of terms of the form [g(x),exp(f(x)),...];
+      -- when 'signum' is +1 (resp. -1), this function filters 'termList'
+      -- leaving only those terms such that coefficient(f(x),ord) is
+      -- maximal (resp. minimal)
+      while (coefficient(exponent first termList,ord) = 0) repeat
+        termList := rest termList
+      empty? termList => error "UPXSSING: can't happen"
+      coefExtreme := coefficient(exponent first termList,ord)
+      outList := list first termList; termList := rest termList
+      for term in termList repeat
+        (coefDiff := coefficient(exponent term,ord) - coefExtreme) = 0 =>
+          outList := concat(term,outList)
+        (sig := sign(coefDiff)$SIGNEF) case "failed" => return "failed"
+        (sig :: Integer) = signum => outList := list term
+      outList
+
+    filterByOrder(termList,predicate) ==
+      -- 'termList' consists of terms of the form [g(x),exp(f(x)),expTerms],
+      -- where 'expTerms' is a list containing some of the terms in the
+      -- series f(x).
+      -- The function filters 'termList' and, when 'predicate' is < (resp. >),
+      -- leaves only those terms with the lowest (resp. highest) order term
+      -- in 'expTerms'
+      while empty? exponentTerms first termList repeat
+        termList := rest termList
+        empty? termList => error "UPXSING: can't happen"
+      ordExtreme := (first exponentTerms first termList).k
+      outList := list first termList
+      for term in rest termList repeat
+        not empty? exponentTerms term =>
+          (ord := (first exponentTerms term).k) = ordExtreme =>
+            outList := concat(term,outList)
+          predicate(ord,ordExtreme) =>
+            ordExtreme := ord
+            outList := list term
+      -- advance pointers on "exponent terms" on terms on 'outList'
+      for term in outList repeat
+        setExponentTerms_!(term,rest exponentTerms term)
+      [outList,ordExtreme]
+
+    dominantTermOnList(termList,ord0,signum) ==
+      -- finds dominant term on 'termList'
+      -- it is known that "exponent terms" of order < 'ord0' are
+      -- the same for all terms on 'termList'
+      newList := termsWithExtremeLeadingCoef(termList,ord0,signum)
+      newList case "failed" => "failed"
+      termList := newList :: List Term
+      empty? rest termList => first termList
+      filtered :=
+        signum = 1 => filterByOrder(termList,#1 < #2)
+        filterByOrder(termList,#1 > #2)
+      termList := filtered.%list
+      empty? rest termList => first termList
+      dominantTermOnList(termList,filtered.%order,signum)
+
+    iDominantTerm termList ==
+      termRecord := sortAndDiscardTerms termList
+      zeroTerms := termRecord.%zeroTerms
+      infiniteTerms := termRecord.%infiniteTerms
+      failedTerms := termRecord.%failedTerms
+      pSeries := termRecord.%puiseuxSeries
+      -- in future versions, we will deal with "failed terms"
+      -- at present, if any occur, we cannot determine the limit
+      not empty? failedTerms => "failed"
+      not zero? pSeries => [makeTerm(pSeries,0),"series"]
+      not empty? infiniteTerms =>
+        empty? rest infiniteTerms => [first infiniteTerms,"infinity"]
+        for term in infiniteTerms repeat computeExponentTerms_! term
+        ord0 := order exponent first infiniteTerms
+        (dTerm := dominantTermOnList(infiniteTerms,ord0,1)) case "failed" =>
+          return "failed"
+        [dTerm :: Term,"infinity"]
+      empty? rest zeroTerms => [first zeroTerms,"zero"]
+      for term in zeroTerms repeat computeExponentTerms_! term
+      ord0 := order exponent first zeroTerms
+      (dTerm := dominantTermOnList(zeroTerms,ord0,-1)) case "failed" =>
+        return "failed"
+      [dTerm :: Term,"zero"]
+
+    dominantTerm f == iDominantTerm terms f
+
+    limitPlus f ==
+      -- list the terms occurring in 'f'; if there are none, then f = 0
+      empty?(termList := terms f) => 0
+      -- compute dominant term
+      (tInfo := iDominantTerm termList) case "failed" => "failed"
+      termInfo := tInfo :: Record(%term:Term,%type:String)
+      domTerm := termInfo.%term
+      (type := termInfo.%type) = "series" =>
+        -- find limit of series term
+        (ord := order(pSeries := coeff domTerm,1)) > 0 => 0
+        coef := coefficient(pSeries,ord)
+        member?(var,variables coef) => "failed"
+        ord = 0 => coef :: OFE
+        -- in the case of an infinite limit, we need to know the sign
+        -- of the first non-zero coefficient
+        (signum := sign(coef)$SIGNEF) case "failed" => "failed"
+        (signum :: Integer) = 1 => plusInfinity()
+        minusInfinity()
+      type = "zero" => 0
+      -- examine lowest order coefficient in series part of 'domTerm'
+      ord := order(pSeries := coeff domTerm)
+      coef := coefficient(pSeries,ord)
+      member?(var,variables coef) => "failed"
+      (signum := sign(coef)$SIGNEF) case "failed" => "failed"
+      (signum :: Integer) = 1 => plusInfinity()
+      minusInfinity()
+
+
+@
+\subsection{gbeuclid.spad}
+\subsubsection{EuclideanGroebnerBasisPackage}
+++ Description: \spadtype{EuclideanGroebnerBasisPackage} computes groebner
+++ bases for polynomial ideals over euclidean domains.
+++ The basic computation provides
+++ a distinguished set of generators for these ideals.
+++ This basis allows an easy test for membership: the operation
+++ \spadfun{euclideanNormalForm} returns zero on ideal members. The string 
+++ "info" and "redcrit" can be given as additional args to provide 
+++ incremental information during the computation. If "info" is given,
+++  a computational summary is given for each s-polynomial. If "redcrit" 
+++ is given, the reduced critical pairs are printed. The term ordering
+++ is determined by the polynomial type used. Suggested types include
+++ \spadtype{DistributedMultivariatePolynomial},
+++ \spadtype{HomogeneousDistributedMultivariatePolynomial},
+++ \spadtype{GeneralDistributedMultivariatePolynomial}.
+<<EuclideanGroebnerBasisPackage>>=
+EuclideanGroebnerBasisPackage(Dom, Expon, VarSet, Dpol): T == C where
+ 
+ Dom: EuclideanDomain
+ Expon: OrderedAbelianMonoidSup
+ VarSet: OrderedSet
+ Dpol: PolynomialCategory(Dom, Expon, VarSet)
+ 
+ T== with
+ 
+     euclideanNormalForm: (Dpol, List(Dpol) )  ->  Dpol
+       ++ euclideanNormalForm(poly,gb) reduces the polynomial poly modulo the
+       ++ precomputed groebner basis gb giving a canonical representative
+       ++ of the residue class.
+     euclideanGroebner: List(Dpol) -> List(Dpol)
+       ++ euclideanGroebner(lp) computes a groebner basis for a polynomial ideal
+       ++ over a euclidean domain generated by the list of polynomials lp.
+     euclideanGroebner: (List(Dpol), String) -> List(Dpol)
+       ++ euclideanGroebner(lp, infoflag) computes a groebner basis 
+       ++ for a polynomial ideal over a euclidean domain
+       ++ generated by the list of polynomials lp.
+       ++ During computation, additional information is printed out
+       ++ if infoflag is given as 
+       ++ either "info" (for summary information) or
+       ++ "redcrit" (for reduced critical pairs)
+     euclideanGroebner: (List(Dpol), String, String ) -> List(Dpol)
+       ++ euclideanGroebner(lp, "info", "redcrit") computes a groebner basis
+       ++ for a polynomial ideal generated by the list of polynomials lp.
+       ++ If the second argument is "info", a summary is given of the critical pairs.
+       ++ If the third argument is "redcrit", critical pairs are printed.
+ C== add
+   Ex ==> OutputForm
+   lc ==> leadingCoefficient
+   red ==> reductum
+
+   import OutputForm
+ 
+   ------  Definition list of critPair
+   ------  lcmfij is now lcm of headterm of poli and polj
+   ------  lcmcij is now lcm of of lc poli and lc polj
+ 
+   critPair ==>Record(lcmfij: Expon, lcmcij: Dom, poli:Dpol, polj: Dpol )
+   Prinp    ==> Record( ci:Dpol,tci:Integer,cj:Dpol,tcj:Integer,c:Dpol,
+                tc:Integer,rc:Dpol,trc:Integer,tH:Integer,tD:Integer)
+ 
+   ------  Definition of intermediate functions
+ 
+   strongGbasis: (List(Dpol), Integer, Integer) -> List(Dpol)
+   eminGbasis: List(Dpol) -> List(Dpol)
+   ecritT: (critPair ) -> Boolean
+   ecritM: (Expon, Dom, Expon, Dom) -> Boolean
+   ecritB: (Expon, Dom, Expon, Dom, Expon, Dom) -> Boolean
+   ecrithinH: (Dpol, List(Dpol)) -> Boolean
+   ecritBonD: (Dpol, List(critPair)) -> List(critPair)
+   ecritMTondd1:(List(critPair)) -> List(critPair)
+   ecritMondd1:(Expon, Dom, List(critPair)) -> List(critPair)
+   crithdelH: (Dpol, List(Dpol)) -> List(Dpol)
+   eupdatF: (Dpol, List(Dpol) ) -> List(Dpol)
+   updatH: (Dpol, List(Dpol), List(Dpol), List(Dpol) ) -> List(Dpol)
+   sortin: (Dpol, List(Dpol) ) -> List(Dpol)
+   eRed: (Dpol, List(Dpol), List(Dpol) )  ->  Dpol
+   ecredPol: (Dpol, List(Dpol) ) -> Dpol
+   esPol: (critPair) -> Dpol
+   updatD: (List(critPair), List(critPair)) -> List(critPair)
+   lepol: Dpol -> Integer
+   prinshINFO : Dpol -> Void
+   prindINFO: (critPair, Dpol, Dpol,Integer,Integer,Integer) -> Integer
+   prinpolINFO: List(Dpol) -> Void
+   prinb: Integer -> Void
+ 
+   ------    MAIN ALGORITHM GROEBNER ------------------------
+   euclideanGroebner( Pol: List(Dpol) ) ==
+     eminGbasis(strongGbasis(Pol,0,0))
+ 
+   euclideanGroebner( Pol: List(Dpol), xx1: String) ==
+     xx1 = "redcrit" =>
+       eminGbasis(strongGbasis(Pol,1,0))
+     xx1 = "info" =>
+       eminGbasis(strongGbasis(Pol,2,1))
+     print("   "::Ex)
+     print("WARNING: options are - redcrit and/or info - "::Ex)
+     print("         you didn't type them correct"::Ex)
+     print("         please try again"::Ex)
+     print("   "::Ex)
+     []
+ 
+   euclideanGroebner( Pol: List(Dpol), xx1: String, xx2: String) ==
+     (xx1 = "redcrit" and xx2 = "info") or
+      (xx1 = "info" and xx2 = "redcrit")   =>
+       eminGbasis(strongGbasis(Pol,1,1))
+     xx1 = "redcrit" and xx2 = "redcrit" =>
+       eminGbasis(strongGbasis(Pol,1,0))
+     xx1 = "info" and xx2 = "info" =>
+       eminGbasis(strongGbasis(Pol,2,1))
+     print("   "::Ex)
+     print("WARNING:  options are - redcrit and/or info - "::Ex)
+     print("          you didn't type them correct"::Ex)
+     print("          please try again "::Ex)
+     print("   "::Ex)
+     []
+ 
+   ------    calculate basis
+ 
+   strongGbasis(Pol: List(Dpol),xx1: Integer, xx2: Integer ) ==
+     dd1, D : List(critPair)
+ 
+     ---------   create D and Pol
+ 
+     Pol1:= sort( (degree #1 > degree #2) or
+                    ((degree #1 = degree #2 ) and
+                       sizeLess?(leadingCoefficient #2,leadingCoefficient #1)),
+                 Pol)
+     Pol:= [first(Pol1)]
+     H:= Pol
+     Pol1:= rest(Pol1)
+     D:= nil
+     while ^null Pol1 repeat
+        h:= first(Pol1)
+        Pol1:= rest(Pol1)
+        en:= degree(h)
+        lch:= lc h
+        dd1:= [[sup(degree(x), en), lcm(leadingCoefficient x, lch), x, h]$critPair
+            for x in Pol]
+        D:= updatD(ecritMTondd1(sort((#1.lcmfij < #2.lcmfij) or
+                                      (( #1.lcmfij = #2.lcmfij ) and
+                                        ( sizeLess?(#1.lcmcij,#2.lcmcij)) ),
+                                     dd1)), ecritBonD(h,D))
+        Pol:= cons(h, eupdatF(h, Pol))
+        ((en = degree(first(H))) and (leadingCoefficient(h) = leadingCoefficient(first(H)) ) ) =>
+              " go to top of while "
+        H:= updatH(h,H,crithdelH(h,H),[h])
+        H:= sort((degree #1 > degree #2) or
+                ((degree #1 = degree #2 ) and
+                  sizeLess?(leadingCoefficient #2,leadingCoefficient #1)), H)
+     D:= sort((#1.lcmfij < #2.lcmfij) or
+             (( #1.lcmfij = #2.lcmfij ) and
+               ( sizeLess?(#1.lcmcij,#2.lcmcij)) ) ,D)
+     xx:= xx2
+ 
+     --------  loop
+ 
+     while ^null D repeat
+         D0:= first D
+         ep:=esPol(D0)
+         D:= rest(D)
+         eh:= ecredPol(eRed(ep,H,H),H)
+         if xx1 = 1 then
+               prinshINFO(eh)
+         eh = 0 =>
+              if xx2 = 1 then
+                  ala:= prindINFO(D0,ep,eh,#H, #D, xx)
+                  xx:= 2
+              " go to top of while "
+         eh := unitCanonical eh
+         e:= degree(eh)
+         leh:= lc eh
+         dd1:= [[sup(degree(x), e), lcm(leadingCoefficient x, leh), x, eh]$critPair
+            for x in Pol]
+         D:= updatD(ecritMTondd1(sort( (#1.lcmfij <
+              #2.lcmfij) or (( #1.lcmfij = #2.lcmfij ) and
+               ( sizeLess?(#1.lcmcij,#2.lcmcij)) ), dd1)), ecritBonD(eh,D))
+         Pol:= cons(eh,eupdatF(eh,Pol))
+         ^ecrithinH(eh,H) or
+           ((e = degree(first(H))) and (leadingCoefficient(eh) = leadingCoefficient(first(H)) ) ) =>
+              if xx2 = 1 then
+                  ala:= prindINFO(D0,ep,eh,#H, #D, xx)
+                  xx:= 2
+              " go to top of while "
+         H:= updatH(eh,H,crithdelH(eh,H),[eh])
+         H:= sort( (degree #1 > degree #2) or
+             ((degree #1 = degree #2 ) and
+                 sizeLess?(leadingCoefficient #2,leadingCoefficient #1)), H)
+         if xx2 = 1 then
+           ala:= prindINFO(D0,ep,eh,#H, #D, xx)
+           xx:= 2
+           " go to top of while "
+     if xx2 = 1 then
+       prinpolINFO(Pol)
+       print("    THE GROEBNER BASIS over EUCLIDEAN DOMAIN"::Ex)
+     if xx1 = 1 and xx2 ^= 1 then
+       print("    THE GROEBNER BASIS over EUCLIDEAN DOMAIN"::Ex)
+     H
+ 
+             --------------------------------------
+ 
+             --- erase multiple of e in D2 using crit M
+ 
+   ecritMondd1(e: Expon, c: Dom, D2: List(critPair))==
+      null D2 => nil
+      x:= first(D2)
+      ecritM(e,c, x.lcmfij, lcm(leadingCoefficient(x.poli), leadingCoefficient(x.polj)))
+         => ecritMondd1(e, c, rest(D2))
+      cons(x, ecritMondd1(e, c, rest(D2)))
+ 
+            -------------------------------
+ 
+   ecredPol(h: Dpol, F: List(Dpol) ) ==
+        h0:Dpol:= 0
+        null F => h
+        while h ^= 0 repeat
+           h0:= h0 + monomial(leadingCoefficient(h),degree(h))
+           h:= eRed(red(h), F, F)
+        h0
+             ----------------------------
+ 
+             --- reduce dd1 using crit T and crit M
+ 
+   ecritMTondd1(dd1: List(critPair))==
+           null dd1 => nil
+           f1:= first(dd1)
+           s1:= #(dd1)
+           cT1:= ecritT(f1)
+           s1= 1 and cT1 => nil
+           s1= 1 => dd1
+           e1:= f1.lcmfij
+           r1:= rest(dd1)
+           f2:= first(r1)
+           e1 = f2.lcmfij and f1.lcmcij = f2.lcmcij =>
+              cT1 =>   ecritMTondd1(cons(f1, rest(r1)))
+              ecritMTondd1(r1)
+           dd1 := ecritMondd1(e1, f1.lcmcij, r1)
+           cT1 => ecritMTondd1(dd1)
+           cons(f1, ecritMTondd1(dd1))
+ 
+             -----------------------------
+ 
+             --- erase elements in D fullfilling crit B
+ 
+   ecritBonD(h:Dpol, D: List(critPair))==
+         null D => nil
+         x:= first(D)
+         x1:= x.poli
+         x2:= x.polj
+         ecritB(degree(h), leadingCoefficient(h), degree(x1),leadingCoefficient(x1),degree(x2),leadingCoefficient(x2)) =>
+           ecritBonD(h, rest(D))
+         cons(x, ecritBonD(h, rest(D)))
+ 
+             -----------------------------
+ 
+             --- concat F and h and erase multiples of h in F
+ 
+   eupdatF(h: Dpol, F: List(Dpol)) ==
+       null F => nil
+       f1:= first(F)
+       ecritM(degree h, leadingCoefficient(h), degree f1, leadingCoefficient(f1))
+           => eupdatF(h, rest(F))
+       cons(f1, eupdatF(h, rest(F)))
+ 
+             -----------------------------
+             --- concat H and h and erase multiples of h in H
+ 
+   updatH(h: Dpol, H: List(Dpol), Hh: List(Dpol), Hhh: List(Dpol)) ==
+       null H => append(Hh,Hhh)
+       h1:= first(H)
+       hlcm:= sup(degree(h1), degree(h))
+       plc:= extendedEuclidean(leadingCoefficient(h), leadingCoefficient(h1))
+       hp:= monomial(plc.coef1,subtractIfCan(hlcm, degree(h))::Expon)*h +
+            monomial(plc.coef2,subtractIfCan(hlcm, degree(h1))::Expon)*h1
+       (ecrithinH(hp, Hh) and ecrithinH(hp, Hhh)) =>
+         hpp:= append(rest(H),Hh)
+         hp:= ecredPol(eRed(hp,hpp,hpp),hpp)
+         updatH(h, rest(H), crithdelH(hp,Hh),cons(hp,crithdelH(hp,Hhh)))
+       updatH(h, rest(H), Hh,Hhh)
+ 
+             --------------------------------------------------
+             ---- delete elements in cons(h,H)
+ 
+   crithdelH(h: Dpol, H: List(Dpol))==
+        null H => nil
+        h1:= first(H)
+        dh1:= degree h1
+        dh:= degree h
+        ecritM(dh, lc h, dh1, lc h1) => crithdelH(h, rest(H))
+        dh1 = sup(dh,dh1) =>
+           plc:= extendedEuclidean( lc h1, lc h)
+           cons(plc.coef1*h1 + monomial(plc.coef2,subtractIfCan(dh1,dh)::Expon)*h,
+               crithdelH(h,rest(H)))
+        cons(h1, crithdelH(h,rest(H)))
+ 
+   eminGbasis(F: List(Dpol)) ==
+        null F => nil
+        newbas := eminGbasis rest F
+        cons(ecredPol( first(F), newbas),newbas)
+ 
+             ------------------------------------------------
+             --- does h belong to H
+ 
+   ecrithinH(h: Dpol, H: List(Dpol))==
+        null H  => true
+        h1:= first(H)
+        ecritM(degree h1, lc h1, degree h, lc h) => false
+        ecrithinH(h, rest(H))
+ 
+            -----------------------------
+            --- calculate  euclidean S-polynomial of a critical pair
+ 
+   esPol(p:critPair)==
+      Tij := p.lcmfij
+      fi := p.poli
+      fj := p.polj
+      lij:= lcm(leadingCoefficient(fi), leadingCoefficient(fj))
+      red(fi)*monomial((lij exquo leadingCoefficient(fi))::Dom,
+                        subtractIfCan(Tij, degree fi)::Expon) -
+        red(fj)*monomial((lij exquo leadingCoefficient(fj))::Dom,
+                         subtractIfCan(Tij, degree fj)::Expon)
+ 
+            ----------------------------
+ 
+            --- euclidean reduction mod F
+ 
+   eRed(s: Dpol, H: List(Dpol), Hh: List(Dpol)) ==
+     ( s = 0 or null H ) => s
+     f1:= first(H)
+     ds:= degree s
+     lf1:= leadingCoefficient(f1)
+     ls:= leadingCoefficient(s)
+     e: Union(Expon, "failed")
+     (((e:= subtractIfCan(ds, degree f1))  case "failed" ) or sizeLess?(ls, lf1) ) =>
+        eRed(s, rest(H), Hh)
+     sdf1:= divide(ls, lf1)
+     q1:= sdf1.quotient
+     sdf1.remainder = 0 =>
+        eRed(red(s) - monomial(q1,e)*reductum(f1), Hh, Hh)
+     eRed(s -(monomial(q1, e)*f1), rest(H), Hh)
+ 
+            ----------------------------
+ 
+            --- crit T  true, if e1 and e2 are disjoint
+ 
+   ecritT(p: critPair) ==
+          pi:= p.poli
+          pj:= p.polj
+          ci:= lc pi
+          cj:= lc pj
+          (p.lcmfij = degree pi + degree pj) and  (p.lcmcij = ci*cj)
+ 
+            ----------------------------
+ 
+            --- crit M - true, if lcm#2 multiple of lcm#1
+ 
+   ecritM(e1: Expon, c1: Dom, e2: Expon, c2: Dom) ==
+     en: Union(Expon, "failed")
+     ((en:=subtractIfCan(e2, e1)) case "failed") or
+       ((c2 exquo c1) case "failed") => false
+     true
+            ----------------------------
+ 
+            --- crit B - true, if eik is a multiple of eh and eik ^equal
+            ---          lcm(eh,ei) and eik ^equal lcm(eh,ek)
+ 
+   ecritB(eh:Expon, ch: Dom, ei:Expon, ci: Dom, ek:Expon, ck: Dom) ==
+       eik:= sup(ei, ek)
+       cik:= lcm(ci, ck)
+       ecritM(eh, ch, eik, cik) and
+             ^ecritM(eik, cik, sup(ei, eh), lcm(ci, ch)) and
+                ^ecritM(eik, cik, sup(ek, eh), lcm(ck, ch))
+ 
+            -------------------------------
+ 
+            --- reduce p1 mod lp
+ 
+   euclideanNormalForm(p1: Dpol, lp: List(Dpol))==
+       eRed(p1, lp, lp)
+ 
+            ---------------------------------
+ 
+            ---  insert element in sorted list
+ 
+   sortin(p1: Dpol, lp: List(Dpol))==
+      null lp => [p1]
+      f1:= first(lp)
+      elf1:= degree(f1)
+      ep1:= degree(p1)
+      ((elf1 < ep1) or ((elf1 = ep1) and
+        sizeLess?(leadingCoefficient(f1),leadingCoefficient(p1)))) =>
+         cons(f1,sortin(p1, rest(lp)))
+      cons(p1,lp)
+ 
+   updatD(D1: List(critPair), D2: List(critPair)) ==
+      null D1 => D2
+      null D2 => D1
+      dl1:= first(D1)
+      dl2:= first(D2)
+      (dl1.lcmfij  <  dl2.lcmfij) => cons(dl1, updatD(D1.rest, D2))
+      cons(dl2, updatD(D1, D2.rest))
+ 
+            ----  calculate number of terms of polynomial
+ 
+   lepol(p1:Dpol)==
+      n: Integer
+      n:= 0
+      while p1 ^= 0 repeat
+         n:= n + 1
+         p1:= red(p1)
+      n
+ 
+            ----  print blanc lines
+ 
+   prinb(n: Integer)==
+        for i in 1..n repeat messagePrint("    ")
+ 
+            ----  print reduced critpair polynom
+ 
+   prinshINFO(h: Dpol)==
+           prinb(2)
+           messagePrint(" reduced Critpair - Polynom :")
+           prinb(2)
+           print(h::Ex)
+           prinb(2)
+ 
+            -------------------------------
+ 
+            ----  print info string
+ 
+   prindINFO(cp: critPair, ps: Dpol, ph: Dpol, i1:Integer,
+             i2:Integer, n:Integer) ==
+       ll: List Prinp
+       a: Dom
+       cpi:= cp.poli
+       cpj:= cp.polj
+       if n = 1 then
+        prinb(1)
+        messagePrint("you choose option  -info-  ")
+        messagePrint("abbrev. for the following information strings are")
+        messagePrint("  ci  =>  Leading monomial  for critpair calculation")
+        messagePrint("  tci =>  Number of terms of polynomial i")
+        messagePrint("  cj  =>  Leading monomial  for critpair calculation")
+        messagePrint("  tcj =>  Number of terms of polynomial j")
+        messagePrint("  c   =>  Leading monomial of critpair polynomial")
+        messagePrint("  tc  =>  Number of terms of critpair polynomial")
+        messagePrint("  rc  =>  Leading monomial of redcritpair polynomial")
+        messagePrint("  trc =>  Number of terms of redcritpair polynomial")
+        messagePrint("  tF  =>  Number of polynomials in reduction list F")
+        messagePrint("  tD  =>  Number of critpairs still to do")
+        prinb(4)
+        n:= 2
+       prinb(1)
+       a:= 1
+       ph = 0  =>
+          ps = 0 =>
+            ll:= [[monomial(a,degree(cpi)),lepol(cpi),monomial(a,degree(cpj)),
+             lepol(cpj),ps,0,ph,0,i1,i2]$Prinp]
+            print(ll::Ex)
+            prinb(1)
+            n
+          ll:= [[monomial(a,degree(cpi)),lepol(cpi),
+            monomial(a,degree(cpj)),lepol(cpj),monomial(a,degree(ps)),
+             lepol(ps), ph,0,i1,i2]$Prinp]
+          print(ll::Ex)
+          prinb(1)
+          n
+       ll:= [[monomial(a,degree(cpi)),lepol(cpi),
+            monomial(a,degree(cpj)),lepol(cpj),monomial(a,degree(ps)),
+             lepol(ps),monomial(a,degree(ph)),lepol(ph),i1,i2]$Prinp]
+       print(ll::Ex)
+       prinb(1)
+       n
+ 
+            -------------------------------
+ 
+            ----  print the groebner basis polynomials
+ 
+   prinpolINFO(pl: List(Dpol))==
+       n:Integer
+       n:= #pl
+       prinb(1)
+       n = 1 =>
+         print("  There is 1  Groebner Basis Polynomial "::Ex)
+         prinb(2)
+       print("  There are "::Ex)
+       prinb(1)
+       print(n::Ex)
+       prinb(1)
+       print("  Groebner Basis Polynomials. "::Ex)
+       prinb(2)
+ 
+@
+\subsection{kl.spad}
+\subsubsection{SortedCache}
+++   A sorted cache of a cachable set S is a dynamic structure that
+++   keeps the elements of S sorted and assigns an integer to each
+++   element of S once it is in the cache. This way, equality and ordering
+++   on S are tested directly on the integers associated with the elements
+++   of S, once they have been entered in the cache.
+<<SortedCache>>=
+SortedCache(S:CachableSet): Exports == Implementation where
+  N    ==> NonNegativeInteger
+  DIFF ==> 1024
+
+  Exports ==> with
+    clearCache  : () -> Void
+      ++ clearCache() empties the cache.
+    cache       : () -> List S
+      ++ cache() returns the current cache as a list.
+    enterInCache: (S, S -> Boolean) -> S
+      ++ enterInCache(x, f) enters x in the cache, calling \spad{f(y)} to
+      ++ determine whether x is equal to y. It returns x with an integer
+      ++ associated with it.
+    enterInCache: (S, (S, S) -> Integer) -> S
+      ++ enterInCache(x, f) enters x in the cache, calling \spad{f(x, y)} to
+      ++ determine whether \spad{x < y (f(x,y) < 0), x = y (f(x,y) = 0)}, or
+      ++ \spad{x > y (f(x,y) > 0)}.
+      ++ It returns x with an integer associated with it.
+
+  Implementation ==> add
+    shiftCache   : (List S, N) -> Void
+    insertInCache: (List S, List S, S, N) -> S
+
+    cach := [nil()]$Record(cche:List S)
+
+    cache() == cach.cche
+
+    shiftCache(l, n) ==
+      for x in l repeat setPosition(x, n + position x)
+      void
+
+    clearCache() ==
+      for x in cache repeat setPosition(x, 0)
+      cach.cche := nil()
+      void
+
+    enterInCache(x:S, equal?:S -> Boolean) ==
+      scan := cache()
+      while not null scan repeat
+        equal?(y := first scan) =>
+          setPosition(x, position y)
+          return y
+        scan := rest scan
+      setPosition(x, 1 + #cache())
+      cach.cche := concat(cache(), x)
+      x
+
+    enterInCache(x:S, triage:(S, S) -> Integer) ==
+      scan := cache()
+      pos:N:= 0
+      for i in 1..#scan repeat
+        zero?(n := triage(x, y := first scan)) =>
+          setPosition(x, position y)
+          return y
+        n<0 => return insertInCache(first(cache(),(i-1)::N),scan,x,pos)
+        scan := rest scan
+        pos  := position y
+      setPosition(x, pos + DIFF)
+      cach.cche := concat(cache(), x)
+      x
+
+    insertInCache(before, after, x, pos) ==
+      if ((pos+1) = position first after) then shiftCache(after, DIFF)
+      setPosition(x, pos + (((position first after) - pos)::N quo 2))
+      cach.cche := concat(before, concat(x, after))
+      x
+
+@
+\subsection{list.spad}
+\subsubsection{IndexedList}
+++   \spadtype{IndexedList} is a basic implementation of the functions
+++   in \spadtype{ListAggregate}, often using functions in the underlying
+++   LISP system. The second parameter to the constructor (\spad{mn})
+++   is the beginning index of the list. That is, if \spad{l} is a
+++   list, then \spad{elt(l,mn)} is the first value. This constructor
+++   is probably best viewed as the implementation of singly-linked
+++   lists that are addressable by index rather than as a mere wrapper
+++   for LISP lists.
+<<IndexedList>>=
+IndexedList(S:Type, mn:Integer): ListAggregate S == add
+   #x                  == LENGTH(x)$Lisp
+   concat(s:S,x:%)     == CONS(s,x)$Lisp
+   eq?(x,y)            == EQ(x,y)$Lisp
+   first x             == SPADfirst(x)$Lisp
+   elt(x,"first")      == SPADfirst(x)$Lisp
+   empty()             == NIL$Lisp
+   empty? x            == NULL(x)$Lisp
+   rest x              == CDR(x)$Lisp
+   elt(x,"rest")       == CDR(x)$Lisp
+   setfirst_!(x,s)     ==
+      empty? x => error "Cannot update an empty list"
+      Qfirst RPLACA(x,s)$Lisp
+   setelt(x,"first",s) ==
+      empty? x => error "Cannot update an empty list"
+      Qfirst RPLACA(x,s)$Lisp
+   setrest_!(x,y)      ==
+      empty? x => error "Cannot update an empty list"
+      Qrest RPLACD(x,y)$Lisp
+   setelt(x,"rest",y)  ==
+      empty? x => error "Cannot update an empty list"
+      Qrest RPLACD(x,y)$Lisp
+   construct l         == l pretend %
+   parts s             == s pretend List S
+   reverse_! x         == NREVERSE(x)$Lisp
+   reverse x           == REVERSE(x)$Lisp
+   minIndex x          == mn
+
+   rest(x, n) ==
+      for i in 1..n repeat
+         if Qnull x then error "index out of range"
+         x := Qrest x
+      x
+
+   copy x ==
+      y := empty()
+      for i in 0.. while not Qnull x repeat
+         if Qeq(i,cycleMax) and cyclic? x then error "cyclic list"
+         y := Qcons(Qfirst x,y)
+         x := Qrest x
+      (NREVERSE(y)$Lisp)@%
+
+   if S has SetCategory then
+     coerce(x):OutputForm ==
+        -- displays cycle with overbar over the cycle
+        y := empty()$List(OutputForm)
+        s := cycleEntry x
+        while Qneq(x, s) repeat
+          y := concat((first x)::OutputForm, y)
+          x := rest x
+        y := reverse_! y
+        empty? s => bracket y
+        -- cyclic case: z is cylic part
+        z := list((first x)::OutputForm)
+        while Qneq(s, rest x) repeat
+           x := rest x
+           z := concat((first x)::OutputForm, z)
+        bracket concat_!(y, overbar commaSeparate reverse_! z)
+
+     x = y ==
+       Qeq(x,y) => true
+       while not Qnull x and not Qnull y repeat
+          Qfirst x ^=$S Qfirst y => return false
+          x := Qrest x
+          y := Qrest y
+       Qnull x and Qnull y
+
+     latex(x : %): String ==
+       s : String := "\left["
+       while not Qnull x repeat
+         s := concat(s, latex(Qfirst x)$S)$String
+         x := Qrest x
+         if not Qnull x then s := concat(s, ", ")$String
+       concat(s, " \right]")$String
+
+     member?(s,x) ==
+        while not Qnull x repeat
+           if s = Qfirst x then return true else x := Qrest x
+        false
+
+   -- Lots of code from parts of AGGCAT, repeated here to
+   -- get faster compilation
+   concat_!(x:%,y:%) ==
+      Qnull x => 
+        Qnull y => x
+        Qpush(first y,x)
+        QRPLACD(x,rest y)$Lisp
+        x
+      z:=x
+      while not Qnull Qrest z repeat
+        z:=Qrest z
+      QRPLACD(z,y)$Lisp
+      x
+
+   -- Then a quicky:
+   if S has SetCategory then
+     removeDuplicates_! l ==
+       p := l
+       while not Qnull p repeat
+--       p := setrest_!(p, remove_!(#1 = Qfirst p, Qrest p))
+-- far too expensive - builds closures etc.
+         pp:=p
+         f:S:=Qfirst p
+         p:=Qrest p
+         while not Qnull (pr:=Qrest pp) repeat
+           if (Qfirst pr)@S = f then QRPLACD(pp,Qrest pr)$Lisp
+           else pp:=pr
+       l
+
+   -- then sorting
+   mergeSort: ((S, S) -> Boolean, %, Integer) -> %
+
+   sort_!(f, l)       == mergeSort(f, l, #l)
+
+   merge_!(f, p, q) ==
+     Qnull p => q
+     Qnull q => p
+     Qeq(p, q) => error "cannot merge a list into itself"
+     if f(Qfirst p, Qfirst q)
+       then (r := t := p; p := Qrest p)
+       else (r := t := q; q := Qrest q)
+     while not Qnull p and not Qnull q repeat
+       if f(Qfirst p, Qfirst q)
+         then (QRPLACD(t, p)$Lisp; t := p; p := Qrest p)
+         else (QRPLACD(t, q)$Lisp; t := q; q := Qrest q)
+     QRPLACD(t, if Qnull p then q else p)$Lisp
+     r
+
+   split_!(p, n) ==
+      n < 1 => error "index out of range"
+      p := rest(p, (n - 1)::NonNegativeInteger)
+      q := Qrest p
+      QRPLACD(p, NIL$Lisp)$Lisp
+      q
+
+   mergeSort(f, p, n) ==
+     if n = 2 and f(first rest p, first p) then p := reverse_! p
+     n < 3 => p
+     l := (n quo 2)::NonNegativeInteger
+     q := split_!(p, l)
+     p := mergeSort(f, p, l)
+     q := mergeSort(f, q, n - l)
+     merge_!(f, p, q)
+
+@
+\subsubsection{List}
+++   \spadtype{List} implements singly-linked lists that are
+++   addressable by indices; the index of the first element
+++   is 1. In addition to the operations provided by
+++   \spadtype{IndexedList}, this constructor provides some
+++   LISP-like functions such as \spadfun{null} and \spadfun{cons}.
+<<List>>=
+List(S:Type): ListAggregate S with
+  nil             : ()     -> %
+    ++ nil() returns the empty list.
+  null            : %      -> Boolean
+    ++ null(u) tests if list \spad{u} is the
+    ++ empty list.
+  cons            : (S, %) -> %
+    ++ cons(element,u) appends \spad{element} onto the front
+    ++ of list \spad{u} and returns the new list. This new list
+    ++ and the old one will share some structure.
+  append          : (%, %) -> %
+    ++ append(u1,u2) appends the elements of list \spad{u1}
+    ++ onto the front of list \spad{u2}. This new list
+    ++ and \spad{u2} will share some structure.
+  if S has SetCategory then
+    setUnion        : (%, %) -> %
+      ++ setUnion(u1,u2) appends the two lists u1 and u2, then
+      ++ removes all duplicates. The order of elements in the
+      ++ resulting list is unspecified.
+    setIntersection : (%, %) -> %
+      ++ setIntersection(u1,u2) returns a list of the elements
+      ++ that lists \spad{u1} and \spad{u2} have in common.
+      ++ The order of elements in the resulting list is unspecified.
+    setDifference   : (%, %) -> %
+      ++ setDifference(u1,u2) returns a list of the elements
+      ++ of \spad{u1} that are not also in \spad{u2}.
+      ++ The order of elements in the resulting list is unspecified.
+  if S has OpenMath then OpenMath
+
+ == IndexedList(S, LISTMININDEX) add
+      nil()                    == NIL$Lisp
+      null l                   == NULL(l)$Lisp
+      cons(s, l)               == CONS(s, l)$Lisp
+      append(l:%, t:%)         == APPEND(l, t)$Lisp
+
+      if S has OpenMath then
+        writeOMList(dev: OpenMathDevice, x: %): Void ==
+          OMputApp(dev)
+          OMputSymbol(dev, "list1", "list")
+          -- The following didn't compile because the compiler isn't
+          -- convinced that `xval' is a S.  Duhhh! MCD.
+          --for xval in x repeat
+          --  OMwrite(dev, xval, false)
+          while not null x repeat
+            OMwrite(dev,first x,false)
+            x := rest x
+          OMputEndApp(dev)
+
+        OMwrite(x: %): String ==
+          s: String := ""
+          sp := OM_-STRINGTOSTRINGPTR(s)$Lisp
+          dev: OpenMathDevice := OMopenString(sp pretend String, OMencodingXML)
+          OMputObject(dev)
+          writeOMList(dev, x)
+          OMputEndObject(dev)
+          OMclose(dev)
+          s := OM_-STRINGPTRTOSTRING(sp)$Lisp pretend String
+          s
+
+        OMwrite(x: %, wholeObj: Boolean): String ==
+          s: String := ""
+          sp := OM_-STRINGTOSTRINGPTR(s)$Lisp
+          dev: OpenMathDevice := OMopenString(sp pretend String, OMencodingXML)
+          if wholeObj then
+            OMputObject(dev)
+          writeOMList(dev, x)
+          if wholeObj then
+            OMputEndObject(dev)
+          OMclose(dev)
+          s := OM_-STRINGPTRTOSTRING(sp)$Lisp pretend String
+          s
+
+        OMwrite(dev: OpenMathDevice, x: %): Void ==
+          OMputObject(dev)
+          writeOMList(dev, x)
+          OMputEndObject(dev)
+
+        OMwrite(dev: OpenMathDevice, x: %, wholeObj: Boolean): Void ==
+          if wholeObj then
+            OMputObject(dev)
+          writeOMList(dev, x)
+          if wholeObj then
+            OMputEndObject(dev)
+
+      if S has SetCategory then
+        setUnion(l1:%,l2:%)      == removeDuplicates concat(l1,l2)
+
+        setIntersection(l1:%,l2:%) ==
+          u :% := empty()
+          l1 := removeDuplicates l1
+          while not empty? l1 repeat
+            if member?(first l1,l2) then u := cons(first l1,u)
+            l1 := rest l1
+          u
+
+        setDifference(l1:%,l2:%) ==
+          l1 := removeDuplicates l1
+          lu:% := empty()
+          while not empty? l1 repeat
+            l11:=l1.1
+            if not member?(l11,l2) then lu := concat(l11,lu)
+            l1 := rest l1
+          lu
+
+      if S has ConvertibleTo InputForm then
+        convert(x:%):InputForm ==
+          convert concat(convert("construct"::Symbol)@InputForm,
+                [convert a for a in (x pretend List S)]$List(InputForm))
+
+
+@
+\subsection{perm.spad}
+\subsubsection{Permutation}
+++ Description: Permutation(S) implements the group of all bijections
+++   on a set S, which move only a finite number of points.
+++   A permutation is considered as a map from S into S. In particular
+++   multiplication is defined as composition of maps:
+++   {\em pi1 * pi2 = pi1 o pi2}.
+++   The internal representation of permuatations are two lists
+++   of equal length representing preimages and images.
+<<Permutation>>=
+Permutation(S:SetCategory): public == private where
+
+  B        ==> Boolean
+  PI       ==> PositiveInteger
+  I        ==> Integer
+  L        ==> List
+  NNI      ==> NonNegativeInteger
+  V        ==> Vector
+  PT       ==> Partition
+  OUTFORM  ==> OutputForm
+  RECCYPE  ==> Record(cycl: L L S, permut: %)
+  RECPRIM  ==> Record(preimage: L S, image : L S)
+
+  public ==> PermutationCategory S with
+
+    listRepresentation:  %                ->  RECPRIM
+      ++ listRepresentation(p) produces a representation {\em rep} of
+      ++ the permutation p as a list of preimages and images, i.e
+      ++ p maps {\em (rep.preimage).k} to {\em (rep.image).k} for all
+      ++ indices k.
+    coercePreimagesImages : List List S   ->  %
+      ++ coercePreimagesImages(lls) coerces the representation {\em lls}
+      ++ of a permutation as a list of preimages and images to a permutation.
+    coerce            :  List List S      ->  %
+      ++ coerce(lls) coerces a list of cycles {\em lls} to a
+      ++ permutation, each cycle being a list with not
+      ++ repetitions, is coerced to the permutation, which maps
+      ++ {\em ls.i} to {\em ls.i+1}, indices modulo the length of the list,
+      ++ then these permutations are mutiplied.
+      ++ Error: if repetitions occur in one cycle.
+    coerce            :  List S           ->  %
+      ++ coerce(ls) coerces a cycle {\em ls}, i.e. a list with not
+      ++ repetitions to a permutation, which maps {\em ls.i} to
+      ++ {\em ls.i+1}, indices modulo the length of the list.
+      ++ Error: if repetitions occur.
+    coerceListOfPairs :  List List S      ->  %
+      ++ coerceListOfPairs(lls) coerces a list of pairs {\em lls} to a
+      ++ permutation.
+      ++ Error: if not consistent, i.e. the set of the first elements
+      ++ coincides with the set of second elements.
+    --coerce            :  %                ->  OUTFORM
+      ++ coerce(p) generates output of the permutation p with domain
+      ++ OutputForm.
+    degree            :  %                ->  NonNegativeInteger
+      ++ degree(p) retuns the number of points moved by the
+      ++ permutation p.
+    movedPoints       :  %                ->  Set S
+      ++ movedPoints(p) returns the set of points moved by the permutation p.
+    cyclePartition    :  %                ->  Partition
+      ++ cyclePartition(p) returns the cycle structure of a permutation
+      ++ p including cycles of length 1 only if S is finite.
+    order             :  %                ->  NonNegativeInteger
+      ++ order(p) returns the order of a permutation p as a group element.
+    numberOfCycles    :  %                ->  NonNegativeInteger
+      ++ numberOfCycles(p) returns the number of non-trivial cycles of
+      ++ the permutation p.
+    sign              :  %                ->  Integer
+      ++ sign(p) returns the signum of the permutation p, +1 or -1.
+    even?             :  %                ->  Boolean
+      ++ even?(p) returns true if and only if p is an even permutation,
+      ++ i.e. {\em sign(p)} is 1.
+    odd?              :  %                ->  Boolean
+      ++ odd?(p) returns true if and only if p is an odd permutation
+      ++ i.e. {\em sign(p)} is {\em -1}.
+    sort              :  L %              ->  L %
+      ++ sort(lp) sorts a list of permutations {\em lp} according to
+      ++ cycle structure first according to length of cycles,
+      ++ second, if S has \spadtype{Finite} or S has
+      ++ \spadtype{OrderedSet} according to lexicographical order of
+      ++ entries in cycles of equal length.
+    if S has Finite then
+      fixedPoints     :  %                ->  Set S
+        ++ fixedPoints(p) returns the points fixed by the permutation p.
+    if S has IntegerNumberSystem or S has Finite then
+      coerceImages    :  L S              ->  %
+        ++ coerceImages(ls) coerces the list {\em ls} to a permutation
+        ++ whose image is given by {\em ls} and the preimage is fixed
+        ++ to be {\em [1,...,n]}.
+        ++ Note: {coerceImages(ls)=coercePreimagesImages([1,...,n],ls)}.
+
+  private ==> add
+
+    -- representation of the object:
+
+    Rep  := V L S
+
+    -- import of domains and packages
+
+    import OutputForm
+    import Vector List S
+
+    -- variables
+
+    p,q      : %
+    exp      : I
+
+    -- local functions first, signatures:
+
+    smaller? : (S,S) -> B
+    rotateCycle: L S -> L S
+    coerceCycle: L L S -> %
+    smallerCycle?: (L S, L S)  -> B
+    shorterCycle?:(L S, L S)  -> B
+    permord:(RECCYPE,RECCYPE)  -> B
+    coerceToCycle:(%,B) -> L L S
+    duplicates?: L S -> B
+
+    smaller?(a:S, b:S): B ==
+      S has OrderedSet => a <$S b
+      S has Finite     => lookup a < lookup b
+      false
+
+    rotateCycle(cyc: L S): L S ==
+      -- smallest element is put in first place
+      -- doesn't change cycle if underlying set
+      -- is not ordered or not finite.
+      min:S := first cyc
+      minpos:I := 1           -- 1 = minIndex cyc
+      for i in 2..maxIndex cyc repeat
+        if smaller?(cyc.i,min) then
+          min  := cyc.i
+          minpos := i
+      one? minpos => cyc
+      concat(last(cyc,((#cyc-minpos+1)::NNI)),first(cyc,(minpos-1)::NNI))
+
+    coerceCycle(lls : L L S): % ==
+      perm : % := 1
+      for lists in reverse lls repeat
+        perm := cycle lists * perm
+      perm
+
+    smallerCycle?(cyca: L S, cycb: L S): B ==
+      #cyca ^= #cycb =>
+        #cyca < #cycb
+      for i in cyca for j in cycb repeat
+        i ^= j => return smaller?(i, j)
+      false
+
+    shorterCycle?(cyca: L S, cycb: L S): B ==
+      #cyca < #cycb
+
+    permord(pa: RECCYPE, pb : RECCYPE): B ==
+      for i in pa.cycl for j in pb.cycl repeat
+        i ^= j => return smallerCycle?(i, j)
+      #pa.cycl < #pb.cycl
+
+    coerceToCycle(p: %, doSorting?: B): L L S ==
+      preim := p.1
+      im := p.2
+      cycles := nil()$(L L S)
+      while not null preim repeat
+        -- start next cycle
+        firstEltInCycle: S := first preim
+        nextCycle : L S := list firstEltInCycle
+        preim := rest preim
+        nextEltInCycle := first im
+        im      := rest im
+        while nextEltInCycle ^= firstEltInCycle repeat
+          nextCycle := cons(nextEltInCycle, nextCycle)
+          i := position(nextEltInCycle, preim)
+          preim := delete(preim,i)
+          nextEltInCycle := im.i
+          im := delete(im,i)
+        nextCycle := reverse nextCycle
+        -- check on 1-cycles, we don't list these
+        if not null rest nextCycle then
+          if doSorting? and (S has OrderedSet or S has Finite) then
+              -- put smallest element in cycle first:
+              nextCycle := rotateCycle nextCycle
+          cycles := cons(nextCycle, cycles)
+      not doSorting? => cycles
+      -- sort cycles
+      S has OrderedSet or S has Finite =>
+        sort(smallerCycle?,cycles)$(L L S)
+      sort(shorterCycle?,cycles)$(L L S)
+
+    duplicates? (ls : L S ): B ==
+      x := copy ls
+      while not null x repeat
+        member? (first x ,rest x) => return true
+        x := rest x
+      false
+
+    -- now the exported functions
+
+    listRepresentation p ==
+      s : RECPRIM := [p.1,p.2]
+
+    coercePreimagesImages preImageAndImage ==
+      p : % := [preImageAndImage.1,preImageAndImage.2]
+
+    movedPoints p == construct p.1  --check on fixed points !!
+
+    degree p ==  #movedPoints p
+
+    p = q ==
+      #(preimp := p.1) ^= #(preimq := q.1) => false
+      for i in 1..maxIndex preimp repeat
+        pos := position(preimp.i, preimq)
+        pos = 0 => return false
+        (p.2).i ^= (q.2).pos => return false
+      true
+
+    orbit(p ,el) ==
+      -- start with a 1-element list:
+      out : Set S := brace list el
+      el2 := eval(p, el)
+      while el2 ^= el repeat
+        -- be carefull: insert adds one element
+        -- as side effect to out
+        insert_!(el2, out)
+        el2 := eval(p, el2)
+      out
+
+    cyclePartition p ==
+      partition([#c for c in coerceToCycle(p, false)])$Partition
+
+    order p ==
+      ord: I := lcm removeDuplicates convert cyclePartition p
+      ord::NNI
+
+    sign(p) ==
+      even? p => 1
+      - 1
+
+
+    even?(p) ==  even?(#(p.1) - numberOfCycles p)
+      -- see the book of James and Kerber on symmetric groups
+      -- for this formula.
+
+    odd?(p) ==  odd?(#(p.1) - numberOfCycles p)
+
+    pa < pb ==
+      pacyc:= coerceToCycle(pa,true)
+      pbcyc:= coerceToCycle(pb,true)
+      for i in pacyc for j in pbcyc repeat
+        i ^= j => return smallerCycle? ( i, j )
+      maxIndex pacyc < maxIndex pbcyc
+
+    coerce(lls : L L S): % == coerceCycle lls
+
+    coerce(ls : L S): % == cycle ls
+
+    sort(inList : L %): L % ==
+      not (S has OrderedSet or S has Finite) => inList
+      ownList: L RECCYPE := nil()$(L RECCYPE)
+      for sigma in inList repeat
+        ownList :=
+          cons([coerceToCycle(sigma,true),sigma]::RECCYPE, ownList)
+      ownList := sort(permord, ownList)$(L RECCYPE)
+      outList := nil()$(L %)
+      for rec in ownList repeat
+        outList := cons(rec.permut, outList)
+      reverse outList
+
+    coerce (p: %): OUTFORM ==
+      cycles: L L S := coerceToCycle(p,true)
+      outfmL : L OUTFORM := nil()
+      for cycle in cycles repeat
+        outcycL: L OUTFORM := nil()
+        for elt in cycle repeat
+          outcycL := cons(elt :: OUTFORM, outcycL)
+        outfmL := cons(paren blankSeparate reverse outcycL, outfmL)
+      -- The identity element will be output as 1:
+      null outfmL => outputForm(1@Integer)
+      -- represent a single cycle in the form (a b c d)
+      -- and not in the form ((a b c d)):
+      null rest outfmL => first outfmL
+      hconcat reverse outfmL
+
+    cycles(vs ) == coerceCycle vs
+
+    cycle(ls) ==
+      #ls < 2 => 1
+      duplicates? ls => error "cycle: the input contains duplicates"
+      [ls, append(rest ls, list first ls)]
+
+    coerceListOfPairs(loP) ==
+      preim := nil()$(L S)
+      im := nil()$(L S)
+      for pair in loP repeat
+        if first pair ^=  second pair then
+          preim := cons(first pair, preim)
+          im := cons(second pair, im)
+      duplicates?(preim) or duplicates?(im) or brace(preim)$(Set S) _
+        ^= brace(im)$(Set S) =>
+        error "coerceListOfPairs: the input cannot be interpreted as a permutation"
+      [preim, im]
+
+    q * p ==
+      -- use vectors for efficiency??
+      preimOfp : V S := construct p.1
+      imOfp : V S := construct p.2
+      preimOfq := q.1
+      imOfq := q.2
+      preimOfqp   := nil()$(L S)
+      imOfqp   := nil()$(L S)
+      -- 1 = minIndex preimOfp
+      for i in 1..(maxIndex preimOfp) repeat
+        -- find index of image of p.i in q if it exists
+        j := position(imOfp.i, preimOfq)
+        if j = 0 then
+          -- it does not exist
+          preimOfqp := cons(preimOfp.i, preimOfqp)
+          imOfqp := cons(imOfp.i, imOfqp)
+        else
+          -- it exists
+          el := imOfq.j
+          -- if the composition fixes the element, we don't
+          -- have to do anything
+          if el ^= preimOfp.i then
+            preimOfqp := cons(preimOfp.i, preimOfqp)
+            imOfqp := cons(el, imOfqp)
+          -- we drop the parts of q which have to do with p
+          preimOfq := delete(preimOfq, j)
+          imOfq := delete(imOfq, j)
+      [append(preimOfqp, preimOfq), append(imOfqp, imOfq)]
+
+    1 == new(2,empty())$Rep
+
+    inv p  == [p.2, p.1]
+
+    eval(p, el) ==
+      pos := position(el, p.1)
+      pos = 0 => el
+      (p.2).pos
+
+    elt(p, el) == eval(p, el)
+
+    numberOfCycles p == #coerceToCycle(p, false)
+
+
+    if S has IntegerNumberSystem then
+
+      coerceImages (image) ==
+        preImage : L S := [i::S for i in 1..maxIndex image]
+        p : % := [preImage,image]
+
+    if S has Finite then
+
+      coerceImages (image) ==
+        preImage : L S := [index(i::PI)::S for i in 1..maxIndex image]
+        p : % := [preImage,image]
+
+      fixedPoints ( p ) == complement movedPoints p
+
+      cyclePartition p ==
+        pt := partition([#c for c in coerceToCycle(p, false)])$Partition
+        pt +$PT conjugate(partition([#fixedPoints(p)])$PT)$PT
+
+@
+\subsection{polset.spad}
+\subsubsection{PolynomialSetCategory}
+<<PolynomialSetCategory>>=
+PolynomialSetCategory(R:Ring, E:OrderedAbelianMonoidSup,_
+  VarSet:OrderedSet, P:RecursivePolynomialCategory(R,E,VarSet)): Category == 
+    Join(SetCategory,Collection(P),CoercibleTo(List(P))) with
+     finiteAggregate
+     retractIfCan : List(P) -> Union($,"failed")
+         ++ \axiom{retractIfCan(lp)} returns an element of the domain whose elements
+         ++ are the members of \axiom{lp} if such an element exists, otherwise
+         ++ \axiom{"failed"} is returned.
+     retract : List(P) -> $
+         ++ \axiom{retract(lp)} returns an element of the domain whose elements
+         ++ are the members of \axiom{lp} if such an element exists, otherwise
+         ++ an error is produced.
+     mvar : $ -> VarSet
+         ++  \axiom{mvar(ps)} returns the main variable of the non constant polynomial
+         ++  with the greatest main variable, if any, else an error is returned.
+     variables : $ -> List VarSet
+         ++  \axiom{variables(ps)} returns the decreasingly sorted list of the
+         ++  variables which are variables of some polynomial in \axiom{ps}.
+     mainVariables : $  -> List VarSet
+         ++ \axiom{mainVariables(ps)} returns the decreasingly sorted list of the
+         ++ variables which are main variables of some polynomial in \axiom{ps}.
+     mainVariable? : (VarSet,$) -> Boolean
+         ++ \axiom{mainVariable?(v,ps)} returns true iff \axiom{v} is the main variable 
+         ++ of some polynomial in \axiom{ps}.
+     collectUnder : ($,VarSet) -> $
+         ++ \axiom{collectUnder(ps,v)} returns the set consisting of the 
+         ++ polynomials of \axiom{ps} with main variable less than \axiom{v}.
+     collect : ($,VarSet) -> $
+         ++ \axiom{collect(ps,v)}  returns the set consisting of the 
+         ++ polynomials of \axiom{ps} with \axiom{v} as main variable.
+     collectUpper : ($,VarSet) -> $
+         ++ \axiom{collectUpper(ps,v)} returns the set consisting of the 
+         ++ polynomials of \axiom{ps} with main variable greater than \axiom{v}.
+     sort : ($,VarSet) -> Record(under:$,floor:$,upper:$)
+         ++ \axiom{sort(v,ps)} returns \axiom{us,vs,ws} such that \axiom{us}
+         ++ is \axiom{collectUnder(ps,v)}, \axiom{vs} is \axiom{collect(ps,v)}
+         ++ and \axiom{ws} is \axiom{collectUpper(ps,v)}. 
+     trivialIdeal?: $ -> Boolean
+         ++ \axiom{trivialIdeal?(ps)} returns true iff \axiom{ps} does
+         ++ not contain non-zero elements.
+     if R has IntegralDomain
+     then
+       roughBase? : $ -> Boolean
+           ++ \axiom{roughBase?(ps)} returns true iff for every pair \axiom{{p,q}}
+           ++ of polynomials in \axiom{ps} their leading monomials are 
+           ++ relatively prime.
+       roughSubIdeal?  : ($,$) -> Boolean
+           ++  \axiom{roughSubIdeal?(ps1,ps2)} returns true iff it can proved 
+           ++ that all polynomials in  \axiom{ps1} lie in the ideal generated by 
+           ++ \axiom{ps2} in  \axiom{\axiom{(R)^(-1) P}} without computing Groebner bases.
+       roughEqualIdeals? : ($,$) -> Boolean
+           ++ \axiom{roughEqualIdeals?(ps1,ps2)} returns true iff it can
+           ++ proved that \axiom{ps1} and \axiom{ps2} generate the same ideal
+           ++ in \axiom{(R)^(-1) P} without computing Groebner bases.
+       roughUnitIdeal? : $ -> Boolean
+           ++ \axiom{roughUnitIdeal?(ps)} returns true iff \axiom{ps} contains some
+           ++ non null element lying in the base ring \axiom{R}.
+       headRemainder : (P,$) -> Record(num:P,den:R)
+           ++ \axiom{headRemainder(a,ps)} returns \axiom{[b,r]} such that the leading
+           ++ monomial of \axiom{b} is reduced in the sense of Groebner bases w.r.t.
+           ++ \axiom{ps} and \axiom{r*a - b} lies in the ideal generated by \axiom{ps}.
+       remainder : (P,$) -> Record(rnum:R,polnum:P,den:R)
+           ++ \axiom{remainder(a,ps)} returns \axiom{[c,b,r]} such that \axiom{b} is fully 
+           ++ reduced in the sense of Groebner bases w.r.t. \axiom{ps}, 
+           ++ \axiom{r*a - c*b} lies in the ideal generated by \axiom{ps}.
+           ++ Furthermore, if \axiom{R} is a gcd-domain, \axiom{b} is primitive.
+       rewriteIdealWithHeadRemainder : (List(P),$) -> List(P)
+           ++ \axiom{rewriteIdealWithHeadRemainder(lp,cs)} returns \axiom{lr} such that
+           ++ the leading monomial of every polynomial in \axiom{lr} is reduced
+           ++ in the sense of Groebner bases w.r.t. \axiom{cs} and \axiom{(lp,cs)}
+           ++ and \axiom{(lr,cs)} generate the same ideal in \axiom{(R)^(-1) P}.
+       rewriteIdealWithRemainder : (List(P),$) -> List(P)
+           ++ \axiom{rewriteIdealWithRemainder(lp,cs)} returns \axiom{lr} such that
+           ++ every polynomial in \axiom{lr} is fully reduced in the sense 
+           ++ of Groebner bases w.r.t. \axiom{cs} and \axiom{(lp,cs)} and 
+           ++ \axiom{(lr,cs)} generate the same ideal in \axiom{(R)^(-1) P}.
+       triangular? : $ -> Boolean
+           ++ \axiom{triangular?(ps)} returns true iff \axiom{ps} is a triangular set,
+           ++ i.e. two distinct polynomials have distinct main variables
+           ++ and no constant lies in \axiom{ps}.
+
+  add
+
+     NNI ==> NonNegativeInteger
+     B ==> Boolean
+
+     elements: $ -> List(P)
+
+     elements(ps:$):List(P) ==
+       lp : List(P) := members(ps)$$
+
+     variables1(lp:List(P)):(List VarSet) ==
+       lvars : List(List(VarSet)) := [variables(p)$P for p in lp]
+       sort(#1 > #2, removeDuplicates(concat(lvars)$List(VarSet)))
+
+     variables2(lp:List(P)):(List VarSet) ==
+       lvars : List(VarSet) := [mvar(p)$P for p in lp]
+       sort(#1 > #2, removeDuplicates(lvars)$List(VarSet))
+
+     variables (ps:$) ==
+       variables1(elements(ps))
+
+     mainVariables (ps:$) ==
+       variables2(remove(ground?,elements(ps)))
+
+     mainVariable? (v,ps) ==
+       lp : List(P) := remove(ground?,elements(ps))
+       while (not empty? lp) and (not (mvar(first(lp)) = v)) repeat
+         lp := rest lp
+       (not empty? lp)
+
+     collectUnder (ps,v) ==
+       lp : List P := elements(ps)
+       lq : List P := []
+       while (not empty? lp) repeat
+         p := first lp
+         lp := rest lp
+         if (ground?(p)) or (mvar(p) < v)
+           then
+             lq := cons(p,lq)
+       construct(lq)$$
+
+     collectUpper (ps,v) ==
+       lp : List P := elements(ps)
+       lq : List P := []
+       while (not empty? lp) repeat
+         p := first lp
+         lp := rest lp
+         if (not ground?(p)) and (mvar(p) > v)
+           then
+             lq := cons(p,lq)
+       construct(lq)$$
+
+     collect (ps,v) ==
+       lp : List P := elements(ps)
+       lq : List P := []
+       while (not empty? lp) repeat
+         p := first lp
+         lp := rest lp
+         if (not ground?(p)) and (mvar(p) = v)
+           then
+             lq := cons(p,lq)
+       construct(lq)$$
+
+     sort (ps,v) ==
+       lp : List P := elements(ps)
+       us : List P := []
+       vs : List P := []
+       ws : List P := []
+       while (not empty? lp) repeat
+         p := first lp
+         lp := rest lp
+         if (ground?(p)) or (mvar(p) < v)
+           then
+             us := cons(p,us)
+           else
+             if (mvar(p) = v)
+               then
+                 vs := cons(p,vs)
+               else
+                 ws := cons(p,ws)
+       [construct(us)$$,construct(vs)$$,construct(ws)$$]$Record(under:$,floor:$,upper:$)
+
+     ps1 = ps2 ==
+       {p for p in elements(ps1)} =$(Set P) {p for p in elements(ps2)}
+
+     exactQuo : (R,R) -> R
+
+     localInf? (p:P,q:P):B ==
+       degree(p) <$E degree(q)
+
+     localTriangular? (lp:List(P)):B ==
+       lp := remove(zero?, lp)
+       empty? lp => true
+       any? (ground?, lp) => false
+       lp := sort(mvar(#1)$P > mvar(#2)$P, lp)
+       p,q : P
+       p := first lp
+       lp := rest lp
+       while (not empty? lp) and (mvar(p) > mvar((q := first(lp)))) repeat
+         p := q
+         lp := rest lp
+       empty? lp
+
+     triangular? ps ==
+       localTriangular? elements ps
+
+     trivialIdeal? ps ==
+       empty?(remove(zero?,elements(ps))$(List(P)))$(List(P))
+
+     if R has IntegralDomain
+     then
+
+       roughUnitIdeal? ps ==
+         any?(ground?,remove(zero?,elements(ps))$(List(P)))$(List P)
+
+       relativelyPrimeLeadingMonomials? (p:P,q:P):B ==
+         dp : E := degree(p)
+         dq : E := degree(q)
+         (sup(dp,dq)$E =$E dp +$E dq)@B
+
+       roughBase? ps ==
+         lp := remove(zero?,elements(ps))$(List(P))
+         empty? lp => true
+         rB? : B := true
+         while (not empty? lp) and rB? repeat
+           p := first lp
+           lp := rest lp
+           copylp := lp
+           while (not empty? copylp) and rB? repeat
+             rB? := relativelyPrimeLeadingMonomials?(p,first(copylp))
+             copylp := rest copylp
+         rB?
+
+       roughSubIdeal?(ps1,ps2) ==
+         lp: List(P) := rewriteIdealWithRemainder(elements(ps1),ps2)
+         empty? (remove(zero?,lp))
+
+       roughEqualIdeals? (ps1,ps2) ==
+         ps1 =$$ ps2 => true
+         roughSubIdeal?(ps1,ps2) and roughSubIdeal?(ps2,ps1)
+
+     if (R has GcdDomain) and (VarSet has ConvertibleTo (Symbol))
+     then
+
+       LPR ==> List Polynomial R
+       LS ==> List Symbol
+
+       if R has EuclideanDomain
+         then
+           exactQuo(r:R,s:R):R ==
+             r quo$R s
+         else
+           exactQuo(r:R,s:R):R ==
+             (r exquo$R s)::R
+
+       headRemainder (a,ps) ==
+         lp1 : List(P) := remove(zero?, elements(ps))$(List(P))
+         empty? lp1 => [a,1$R]
+         any?(ground?,lp1) => [reductum(a),1$R]
+         r : R := 1$R
+         lp1 := sort(localInf?, reverse elements(ps))
+         lp2 := lp1
+         e : Union(E, "failed")
+         while (not zero? a) and (not empty? lp2) repeat
+           p := first lp2
+           if ((e:= subtractIfCan(degree(a),degree(p))) case E)
+             then
+               g := gcd((lca := leadingCoefficient(a)),(lcp := leadingCoefficient(p)))$R
+               (lca,lcp) := (exactQuo(lca,g),exactQuo(lcp,g))
+               a := lcp * reductum(a) - monomial(lca, e::E)$P * reductum(p)
+               r := r * lcp
+               lp2 := lp1
+             else
+               lp2 := rest lp2
+         [a,r]
+
+       makeIrreducible! (frac:Record(num:P,den:R)):Record(num:P,den:R) ==
+         g := gcd(frac.den,frac.num)$P
+         one? g => frac
+         frac.num := exactQuotient!(frac.num,g)
+         frac.den := exactQuo(frac.den,g)
+         frac
+
+       remainder (a,ps) ==
+         hRa := makeIrreducible! headRemainder (a,ps)
+         a := hRa.num
+         r : R := hRa.den
+         zero? a => [1$R,a,r]
+         b : P := monomial(1$R,degree(a))$P
+         c : R := leadingCoefficient(a)
+         while not zero?(a := reductum a) repeat
+           hRa := makeIrreducible!  headRemainder (a,ps)
+           a := hRa.num
+           r := r * hRa.den
+           g := gcd(c,(lca := leadingCoefficient(a)))$R
+           b := ((hRa.den) * exactQuo(c,g)) * b + monomial(exactQuo(lca,g),degree(a))$P
+           c := g
+         [c,b,r]
+
+       rewriteIdealWithHeadRemainder(ps,cs) ==
+         trivialIdeal? cs => ps
+         roughUnitIdeal? cs => [0$P]
+         ps := remove(zero?,ps)
+         empty? ps => ps
+         any?(ground?,ps) => [1$P]
+         rs : List P := []
+         while not empty? ps repeat
+           p := first ps
+           ps := rest ps
+           p := (headRemainder(p,cs)).num
+           if not zero? p
+             then 
+               if ground? p
+                 then
+                   ps := []
+                   rs := [1$P]
+                 else
+                   primitivePart! p
+                   rs := cons(p,rs)
+         removeDuplicates rs
+
+       rewriteIdealWithRemainder(ps,cs) ==
+         trivialIdeal? cs => ps
+         roughUnitIdeal? cs => [0$P]
+         ps := remove(zero?,ps)
+         empty? ps => ps
+         any?(ground?,ps) => [1$P]
+         rs : List P := []
+         while not empty? ps repeat
+           p := first ps
+           ps := rest ps
+           p := (remainder(p,cs)).polnum
+           if not zero? p
+             then 
+               if ground? p
+                 then
+                   ps := []
+                   rs := [1$P]
+                 else
+                   rs := cons(unitCanonical(p),rs)
+         removeDuplicates rs
+
+@
+\subsection{sortpak.spad}
+\subsubsection{SortPackage}
+++ This package exports sorting algorithnms
+<<SortPackage>>=
+SortPackage(S,A) : Exports == Implementation where
+  S: Type
+  A: IndexedAggregate(Integer,S)
+    with (finiteAggregate; shallowlyMutable)
+
+  Exports == with
+    bubbleSort_!: (A,(S,S) -> Boolean) -> A
+	++ bubbleSort!(a,f) \undocumented
+    insertionSort_!: (A, (S,S) -> Boolean) -> A
+	++ insertionSort!(a,f) \undocumented
+    if S has OrderedSet then
+      bubbleSort_!: A -> A
+	++ bubbleSort!(a) \undocumented
+      insertionSort_!: A -> A
+	++ insertionSort! \undocumented
+
+  Implementation == add
+    bubbleSort_!(m,f) ==
+      n := #m
+      for i in 1..(n-1) repeat
+        for j in n..(i+1) by -1 repeat
+          if f(m.j,m.(j-1)) then swap_!(m,j,j-1)
+      m
+    insertionSort_!(m,f) ==
+      for i in 2..#m repeat
+        j := i
+        while j > 1 and f(m.j,m.(j-1)) repeat
+          swap_!(m,j,j-1)
+          j := (j - 1) pretend PositiveInteger
+      m
+    if S has OrderedSet then
+      bubbleSort_!(m) == bubbleSort_!(m,_<$S)
+      insertionSort_!(m) == insertionSort_!(m,_<$S)
+    if A has UnaryRecursiveAggregate(S) then
+      bubbleSort_!(m,fn) ==
+        empty? m => m
+        l := m
+        while not empty? (r := l.rest) repeat
+           r := bubbleSort_!(r,fn)
+           x := l.first
+           if fn(r.first,x) then
+             l.first := r.first
+             r.first := x
+           l.rest := r
+           l := l.rest
+        m
+
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} {\bf PseudoRemainderSequence in prs.spad}
+\end{thebibliography}
+\end{document}