diff options
Diffstat (limited to 'src/algebra/lodof.spad.pamphlet')
| -rw-r--r-- | src/algebra/lodof.spad.pamphlet | 533 |
1 files changed, 533 insertions, 0 deletions
diff --git a/src/algebra/lodof.spad.pamphlet b/src/algebra/lodof.spad.pamphlet new file mode 100644 index 00000000..fd72c15a --- /dev/null +++ b/src/algebra/lodof.spad.pamphlet @@ -0,0 +1,533 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra lodof.spad} +\author{Manuel Bronstein, Fritz Schwarz} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{domain SETMN SetOfMIntegersInOneToN} +<<domain SETMN SetOfMIntegersInOneToN>>= +)abbrev domain SETMN SetOfMIntegersInOneToN +++ Author: Manuel Bronstein +++ Date Created: 10 January 1994 +++ Date Last Updated: 10 January 1994 +++ Description: +++ \spadtype{SetOfMIntegersInOneToN} implements the subsets of M integers +++ in the interval \spad{[1..n]} +SetOfMIntegersInOneToN(m, n): Exports == Implementation where + PI ==> PositiveInteger + N ==> NonNegativeInteger + U ==> Union(%, "failed") + n,m: PI + + Exports ==> Finite with + incrementKthElement: (%, PI) -> U + ++ incrementKthElement(S,k) increments the k^{th} element of S, + ++ and returns "failed" if the result is not a set of M integers + ++ in \spad{1..n} any more. + replaceKthElement: (%, PI, PI) -> U + ++ replaceKthElement(S,k,p) replaces the k^{th} element of S by p, + ++ and returns "failed" if the result is not a set of M integers + ++ in \spad{1..n} any more. + elements: % -> List PI + ++ elements(S) returns the list of the elements of S in increasing order. + setOfMinN: List PI -> % + ++ setOfMinN([a_1,...,a_m]) returns the set {a_1,...,a_m}. + ++ Error if {a_1,...,a_m} is not a set of M integers in \spad{1..n}. + enumerate: () -> Vector % + ++ enumerate() returns a vector of all the sets of M integers in + ++ \spad{1..n}. + member?: (PI, %) -> Boolean + ++ member?(p, s) returns true is p is in s, false otherwise. + delta: (%, PI, PI) -> N + ++ delta(S,k,p) returns the number of elements of S which are strictly + ++ between p and the k^{th} element of S. + + Implementation ==> add + Rep := Record(bits:Bits, pos:N) + + reallyEnumerate: () -> Vector % + enum: (N, N, PI) -> List Bits + + all:Reference Vector % := ref empty() + sz:Reference N := ref 0 + + s1 = s2 == s1.bits =$Bits s2.bits + coerce(s:%):OutputForm == brace [i::OutputForm for i in elements s] + random() == index((1 + (random()$Integer rem size()))::PI) + reallyEnumerate() == [[b, i] for b in enum(m, n, n) for i in 1..] + member?(p, s) == s.bits.p + + enumerate() == + if empty? all() then all() := reallyEnumerate() + all() + +-- enumerates the sets of p integers in 1..q, returns them as sets in 1..n +-- must have p <= q + enum(p, q, n) == + zero? p or zero? q => empty() + p = q => + b := new(n, false)$Bits + for i in 1..p repeat b.i := true + [b] + q1 := (q - 1)::N + l := enum((p - 1)::N, q1, n) + if empty? l then l := [new(n, false)$Bits] + for s in l repeat s.q := true + concat_!(enum(p, q1, n), l) + + size() == + if zero? sz() then + sz() := binomial(n, m)$IntegerCombinatoricFunctions(Integer) :: N + sz() + + lookup s == + if empty? all() then all() := reallyEnumerate() + if zero?(s.pos) then s.pos := position(s, all()) :: N + s.pos :: PI + + index p == + p > size() => error "index: argument too large" + if empty? all() then all() := reallyEnumerate() + all().p + + setOfMinN l == + s := new(n, false)$Bits + count:N := 0 + for i in l repeat + count := count + 1 + count > m or zero? i or i > n or s.i => + error "setOfMinN: improper set of integers" + s.i := true + count < m => error "setOfMinN: improper set of integers" + [s, 0] + + elements s == + b := s.bits + l:List PI := empty() + found:N := 0 + i:PI := 1 + while found < m repeat + if b.i then + l := concat(i, l) + found := found + 1 + i := i + 1 + reverse_! l + + incrementKthElement(s, k) == + b := s.bits + found:N := 0 + i:N := 1 + while found < k repeat + if b.i then found := found + 1 + i := i + 1 + i > n or b.i => "failed" + newb := copy b + newb.i := true + newb.((i-1)::N) := false + [newb, 0] + + delta(s, k, p) == + b := s.bits + count:N := found:N := 0 + i:PI := 1 + while found < k repeat + if b.i then + found := found + 1 + if i > p and found < k then count := count + 1 + i := i + 1 + count + + replaceKthElement(s, k, p) == + b := s.bits + found:N := 0 + i:PI := 1 + while found < k repeat + if b.i then found := found + 1 + if found < k then i := i + 1 + b.p and i ^= p => "failed" + newb := copy b + newb.p := true + newb.i := false + [newb, (i = p => s.pos; 0)] + +@ +\section{package PREASSOC PrecomputedAssociatedEquations} +<<package PREASSOC PrecomputedAssociatedEquations>>= +)abbrev package PREASSOC PrecomputedAssociatedEquations +++ Author: Manuel Bronstein +++ Date Created: 13 January 1994 +++ Date Last Updated: 3 February 1994 +++ Description: +++ \spadtype{PrecomputedAssociatedEquations} stores some generic +++ precomputations which speed up the computations of the +++ associated equations needed for factoring operators. +PrecomputedAssociatedEquations(R, L): Exports == Implementation where + R: IntegralDomain + L: LinearOrdinaryDifferentialOperatorCategory R + + PI ==> PositiveInteger + N ==> NonNegativeInteger + A ==> PrimitiveArray R + U ==> Union(Matrix R, "failed") + + Exports ==> with + firstUncouplingMatrix: (L, PI) -> U + ++ firstUncouplingMatrix(op, m) returns the matrix A such that + ++ \spad{A w = (W',W'',...,W^N)} in the corresponding associated + ++ equations for right-factors of order m of op. + ++ Returns "failed" if the matrix A has not been precomputed for + ++ the particular combination \spad{degree(L), m}. + + Implementation ==> add + A32: L -> U + A42: L -> U + A425: (A, A, A) -> List R + A426: (A, A, A) -> List R + makeMonic: L -> Union(A, "failed") + + diff:L := D() + + firstUncouplingMatrix(op, m) == + n := degree op + n = 3 and m = 2 => A32 op + n = 4 and m = 2 => A42 op + "failed" + + makeMonic op == + lc := leadingCoefficient op + a:A := new(n := degree op, 0) + for i in 0..(n-1)::N repeat + (u := coefficient(op, i) exquo lc) case "failed" => return "failed" + a.i := - (u::R) + a + + A32 op == + (u := makeMonic op) case "failed" => "failed" + a := u::A + matrix [[0, 1, 0], [a.1, a.2, 1], + [diff(a.1) + a.1 * a.2 - a.0, diff(a.2) + a.2**2 + a.1, 2 * a.2]] + + A42 op == + (u := makeMonic op) case "failed" => "failed" + a := u::A + a':A := new(4, 0) + a'':A := new(4, 0) + for i in 0..3 repeat + a'.i := diff(a.i) + a''.i := diff(a'.i) + matrix [[0, 1, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0], [a.1,a.2,0,a.3,2::R,0], + [a'.1 + a.1 * a.3 - 2 * a.0, a'.2 + a.2 * a.3 + a.1, 3 * a.2, + a'.3 + a.3 ** 2 + a.2, 3 * a.3, 2::R], + A425(a, a', a''), A426(a, a', a'')] + + A425(a, a', a'') == + [a''.1 + 2 * a.1 * a'.3 + a.3 * a'.1 - 2 * a'.0 + a.1 * a.3 ** 2 + - 3 * a.0 * a.3 + a.1 * a.2, + a''.2 + 2 * a.2 * a'.3 + a.3 * a'.2 + 2 * a'.1 + a.2 * a.3 ** 2 + + a.1 * a.3 + a.2 ** 2 - 4 * a.0, + 4 * a'.2 + 4 * a.2 * a.3 - a.1, + a''.3 + 3 * a.3 * a'.3 + 2 * a'.2 + a.3 ** 3 + 2 * a.2 * a.3 + a.1, + 4 * a'.3 + 4 * a.3 ** 2 + 4 * a.2, 5 * a.3] + + A426(a, a', a'') == + [diff(a''.1) + 3 * a.1 * a''.3 + a.3 * a''.1 - 2 * a''.0 + + (3 * a'.1 + 5 * a.1 * a.3 - 7 * a.0) * a'.3 + 3 * a.1 * a'.2 + + (a.3 ** 2 + a.2) * a'.1 - 3 * a.3 * a'.0 + a.1 * a.3 ** 3 + - 4 * a.0 * a.3 ** 2 + 2 * a.1 * a.2 * a.3 - 4 * a.0 * a.2 + a.1 ** 2, + diff(a''.2) + 3 * a.2 * a''.3 + a.3 * a''.2 + 3 * a''.1 + + (3*a'.2 + 5*a.2 * a.3 + 3 * a.1) * a'.3 + (a.3**2 + 4*a.2)*a'.2 + + 2 * a.3 * a'.1 - 6 * a'.0 + a.2 * a.3 ** 3 + a.1 * a.3 ** 2 + + (2 * a.2**2 - 8 * a.0) * a.3 + 2 * a.1 * a.2, + 5 * a''.2 + 10 * a.2 * a'.3 + 5 * a.3 * a'.2 + a'.1 + + 5 * a.2 * a.3 ** 2 - 4 * a.1 * a.3 + 5 * a.2**2 - 4 * a.0, + diff(a''.3) + 4 * a.3 * a''.3 + 3*a''.2 + 3 * a'.3**2 + + (6 * a.3**2 + 4 * a.2) * a'.3 + 5 * a.3 * a'.2 + 3 * a'.1 + + a.3**4 + 3 * a.2 * a.3**2 + 2 * a.1 * a.3 + a.2**2 - 4*a.0, + 5 * a''.3 + 15 * a.3 * a'.3 + 10 * a'.2 + 5 * a.3**3 + + 10 * a.2 * a.3, 9 * a'.3 + 9 * a.3**2 + 4 * a.2] + +@ +\section{package ASSOCEQ AssociatedEquations} +<<package ASSOCEQ AssociatedEquations>>= +)abbrev package ASSOCEQ AssociatedEquations +++ Author: Manuel Bronstein +++ Date Created: 10 January 1994 +++ Date Last Updated: 3 February 1994 +++ Description: +++ \spadtype{AssociatedEquations} provides functions to compute the +++ associated equations needed for factoring operators +AssociatedEquations(R, L):Exports == Implementation where + R: IntegralDomain + L: LinearOrdinaryDifferentialOperatorCategory R + + PI ==> PositiveInteger + N ==> NonNegativeInteger + MAT ==> Matrix R + REC ==> Record(minor: List PI, eq: L, minors: List List PI, ops: List L) + + Exports ==> with + associatedSystem: (L, PI) -> Record(mat: MAT, vec:Vector List PI) + ++ associatedSystem(op, m) returns \spad{[M,w]} such that the + ++ m-th associated equation system to L is \spad{w' = M w}. + uncouplingMatrices: MAT -> Vector MAT + ++ uncouplingMatrices(M) returns \spad{[A_1,...,A_n]} such that if + ++ \spad{y = [y_1,...,y_n]} is a solution of \spad{y' = M y}, then + ++ \spad{[$y_j',y_j'',...,y_j^{(n)}$] = $A_j y$} for all j's. + if R has Field then + associatedEquations: (L, PI) -> REC + ++ associatedEquations(op, m) returns \spad{[w, eq, lw, lop]} + ++ such that \spad{eq(w) = 0} where w is the given minor, and + ++ \spad{lw_i = lop_i(w)} for all the other minors. + + Implementation ==> add + makeMatrix: (Vector MAT, N) -> MAT + + diff:L := D() + + makeMatrix(v, n) == matrix [parts row(v.i, n) for i in 1..#v] + + associatedSystem(op, m) == + eq: Vector R + S := SetOfMIntegersInOneToN(m, n := degree(op)::PI) + w := enumerate()$S + s := size()$S + ww:Vector List PI := new(s, empty()) + M:MAT := new(s, s, 0) + m1 := (m::Integer - 1)::PI + an := leadingCoefficient op + a:Vector(R) := [- (coefficient(op, j) exquo an)::R for j in 0..n - 1] + for i in 1..s repeat + eq := new(s, 0) + wi := w.i + ww.i := elements wi + for k in 1..m1 repeat + u := incrementKthElement(wi, k::PI)$S + if u case S then eq(lookup(u::S)) := 1 + if member?(n, wi) then + for j in 1..n | a.j ^= 0 repeat + u := replaceKthElement(wi, m, j::PI) + if u case S then + eq(lookup(u::S)) := (odd? delta(wi, m, j::PI) => -a.j; a.j) + else + u := incrementKthElement(wi, m)$S + if u case S then eq(lookup(u::S)) := 1 + setRow_!(M, i, eq) + [M, ww] + + uncouplingMatrices m == + n := nrows m + v:Vector MAT := new(n, zero(1, 0)$MAT) + v.1 := mi := m + for i in 2..n repeat v.i := mi := map(diff #1, mi) + mi * m + [makeMatrix(v, i) for i in 1..n] + + if R has Field then + import PrecomputedAssociatedEquations(R, L) + + makeop: Vector R -> L + makeeq: (Vector List PI, MAT, N, N) -> REC + computeIt: (L, PI, N) -> REC + + makeeq(v, m, i, n) == + [v.i, makeop row(m, i) - 1, [v.j for j in 1..n | j ^= i], + [makeop row(m, j) for j in 1..n | j ^= i]] + + associatedEquations(op, m) == + (u := firstUncouplingMatrix(op, m)) case "failed" => computeIt(op,m,1) + (v := inverse(u::MAT)) case "failed" => computeIt(op, m, 2) + S := SetOfMIntegersInOneToN(m, degree(op)::PI) + w := enumerate()$S + s := size()$S + ww:Vector List PI := new(s, empty()) + for i in 1..s repeat ww.i := elements(w.i) + makeeq(ww, v::MAT, 1, s) + + computeIt(op, m, k) == + rec := associatedSystem(op, m) + a := uncouplingMatrices(rec.mat) + n := #a + for i in k..n repeat + (u := inverse(a.i)) case MAT => return makeeq(rec.vec,u::MAT,i,n) + error "associatedEquations: full degenerate case" + + makeop v == + op:L := 0 + for i in 1..#v repeat op := op + monomial(v i, i) + op + +@ +\section{package LODOF LinearOrdinaryDifferentialOperatorFactorizer} +<<package LODOF LinearOrdinaryDifferentialOperatorFactorizer>>= +)abbrev package LODOF LinearOrdinaryDifferentialOperatorFactorizer +++ Author: Fritz Schwarz, Manuel Bronstein +++ Date Created: 1988 +++ Date Last Updated: 3 February 1994 +++ Description: +++ \spadtype{LinearOrdinaryDifferentialOperatorFactorizer} provides a +++ factorizer for linear ordinary differential operators whose coefficients +++ are rational functions. +++ Keywords: differential equation, ODE, LODO, factoring +LinearOrdinaryDifferentialOperatorFactorizer(F, UP): Exports == Impl where + F : Join(Field, CharacteristicZero, + RetractableTo Integer, RetractableTo Fraction Integer) + UP: UnivariatePolynomialCategory F + + RF ==> Fraction UP + L ==> LinearOrdinaryDifferentialOperator1 RF + + Exports ==> with + factor: (L, UP -> List F) -> List L + ++ factor(a, zeros) returns the factorisation of a. + ++ \spad{zeros} is a zero finder in \spad{UP}. + if F has AlgebraicallyClosedField then + factor: L -> List L + ++ factor(a) returns the factorisation of a. + factor1: L -> List L + ++ factor1(a) returns the factorisation of a, + ++ assuming that a has no first-order right factor. + + Impl ==> add + import RationalLODE(F, UP) + import RationalRicDE(F, UP) +-- import AssociatedEquations RF + + dd := D()$L + + expsol : (L, UP -> List F, UP -> Factored UP) -> Union(RF, "failed") + expsols : (L, UP -> List F, UP -> Factored UP, Boolean) -> List RF + opeval : (L, L) -> L + recurfactor: (L, L, UP -> List F, UP -> Factored UP, Boolean) -> List L + rfactor : (L, L, UP -> List F, UP -> Factored UP, Boolean) -> List L + rightFactor: (L, NonNegativeInteger, UP -> List F, UP -> Factored UP) + -> Union(L, "failed") + innerFactor: (L, UP -> List F, UP -> Factored UP, Boolean) -> List L + + factor(l, zeros) == innerFactor(l, zeros, squareFree, true) + + expsol(l, zeros, ezfactor) == + empty?(sol := expsols(l, zeros, ezfactor, false)) => "failed" + first sol + + expsols(l, zeros, ezfactor, all?) == + sol := [differentiate(f)/f for f in ratDsolve(l, 0).basis | f ^= 0] + not(all? or empty? sol) => sol + concat(sol, ricDsolve(l, zeros, ezfactor)) + +-- opeval(l1, l2) returns l1(l2) + opeval(l1, l2) == + ans:L := 0 + l2n:L := 1 + for i in 0..degree l1 repeat + ans := ans + coefficient(l1, i) * l2n + l2n := l2 * l2n + ans + + recurfactor(l, r, zeros, ezfactor, adj?) == + q := rightExactQuotient(l, r)::L + if adj? then q := adjoint q + innerFactor(q, zeros, ezfactor, true) + + rfactor(op, r, zeros, ezfactor, adj?) == +-- degree r > 1 or not one? leadingCoefficient r => + degree r > 1 or not ((leadingCoefficient r) = 1) => + recurfactor(op, r, zeros, ezfactor, adj?) + op1 := opeval(op, dd - coefficient(r, 0)::L) + map_!(opeval(#1, r), recurfactor(op1, dd, zeros, ezfactor, adj?)) + +-- r1? is true means look for 1st-order right-factor also + innerFactor(l, zeros, ezfactor, r1?) == + (n := degree l) <= 1 => [l] + ll := adjoint l + for i in 1..(n quo 2) repeat + (r1? or (i > 1)) and ((u := rightFactor(l,i,zeros,ezfactor)) case L) => + return concat_!(rfactor(l, u::L, zeros, ezfactor, false), u::L) + (2 * i < n) and ((u := rightFactor(ll, i, zeros, ezfactor)) case L) => + return concat(adjoint(u::L), rfactor(ll, u::L, zeros,ezfactor,true)) + [l] + + rightFactor(l, n, zeros, ezfactor) == +-- one? n => + (n = 1) => + (u := expsol(l, zeros, ezfactor)) case "failed" => "failed" + D() - u::RF::L +-- rec := associatedEquations(l, n::PositiveInteger) +-- empty?(sol := expsols(rec.eq, zeros, ezfactor, true)) => "failed" + "failed" + + if F has AlgebraicallyClosedField then + zro1: UP -> List F + zro : (UP, UP -> Factored UP) -> List F + + zro(p, ezfactor) == + concat [zro1(r.factor) for r in factors ezfactor p] + + zro1 p == + [zeroOf(map(#1, p)$UnivariatePolynomialCategoryFunctions2(F, UP, + F, SparseUnivariatePolynomial F))] + + if F is AlgebraicNumber then + import AlgFactor UP + + factor l == innerFactor(l, zro(#1, factor), factor, true) + factor1 l == innerFactor(l, zro(#1, factor), factor, false) + + else + factor l == innerFactor(l, zro(#1, squareFree), squareFree, true) + factor1 l == innerFactor(l, zro(#1, squareFree), squareFree, false) + +@ +\section{License} +<<license>>= +--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. +--All rights reserved. +-- +--Redistribution and use in source and binary forms, with or without +--modification, are permitted provided that the following conditions are +--met: +-- +-- - Redistributions of source code must retain the above copyright +-- notice, this list of conditions and the following disclaimer. +-- +-- - Redistributions in binary form must reproduce the above copyright +-- notice, this list of conditions and the following disclaimer in +-- the documentation and/or other materials provided with the +-- distribution. +-- +-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the +-- names of its contributors may be used to endorse or promote products +-- derived from this software without specific prior written permission. +-- +--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +@ +<<*>>= +<<license>> + +-- Compile order for the differential equation solver: +-- oderf.spad odealg.spad nlode.spad nlinsol.spad riccati.spad +-- kovacic.spad lodof.spad odeef.spad + +<<domain SETMN SetOfMIntegersInOneToN>> +<<package PREASSOC PrecomputedAssociatedEquations>> +<<package ASSOCEQ AssociatedEquations>> +<<package LODOF LinearOrdinaryDifferentialOperatorFactorizer>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |