\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(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(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(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(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(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(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(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(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(): 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(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(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(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(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(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(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(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(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}