diff options
| -rw-r--r-- | src/ChangeLog | 4 | ||||
| -rw-r--r-- | src/booklets/ChangeLog | 18 | ||||
| -rw-r--r-- | src/booklets/Makefile.in | 15 | ||||
| -rw-r--r-- | src/booklets/Makefile.pamphlet | 36 | ||||
| -rwxr-xr-x | src/booklets/Rosetta.booklet | 90 | ||||
| -rwxr-xr-x | src/booklets/Sorting.booklet | 2731 |
6 files changed, 4 insertions, 2890 deletions
diff --git a/src/ChangeLog b/src/ChangeLog index bbd67531..3819c0f3 100644 --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,7 @@ +2007-12-27 Gabriel Dos Reis <gdr@cs.tamu.edu> + + * booklets: Remove. + 2007-11-25 Gabriel Dos Reis <gdr@cs.tamu.edu> * Makefile.pamphlet: Remove all-depsys rule. diff --git a/src/booklets/ChangeLog b/src/booklets/ChangeLog deleted file mode 100644 index 305e6d9b..00000000 --- a/src/booklets/ChangeLog +++ /dev/null @@ -1,18 +0,0 @@ -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 deleted file mode 100644 index 97aae83e..00000000 --- a/src/booklets/Makefile.in +++ /dev/null @@ -1,15 +0,0 @@ - -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 deleted file mode 100644 index d0573f47..00000000 --- a/src/booklets/Makefile.pamphlet +++ /dev/null @@ -1,36 +0,0 @@ -%% 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 deleted file mode 100755 index 171c0f0e..00000000 --- a/src/booklets/Rosetta.booklet +++ /dev/null @@ -1,90 +0,0 @@ -\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 deleted file mode 100755 index 2f43240f..00000000 --- a/src/booklets/Sorting.booklet +++ /dev/null @@ -1,2731 +0,0 @@ -\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} |