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