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