path: root/src/booklets/Sorting.booklet

\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}