\documentclass{article} \usepackage{open-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} <>= )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) -> 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(size())$Integer)::PI) reallyEnumerate() == [[b, i] for b in enum(m, n) for i in 1..] member?(p, s) == s.bits.p enumerate():Vector % == if empty? deref all then setref(all,reallyEnumerate()) deref all -- enumerates the sets of p integers in 1..q, returns them as sets in 1..n -- must have p <= q enum(p, q) == 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) if empty? l then l := [new(n, false)$Bits] for s in l repeat s.q := true concat!(enum(p, q1), l) size() == if zero? deref sz then setref(sz,binomial(n, m)$IntegerCombinatoricFunctions(Integer) :: N) deref sz lookup s == if empty? deref all then setref(all,reallyEnumerate()) if zero?(s.pos) then s.pos := position(s, deref all) :: N s.pos :: PI index p == p > size() => error "index: argument too large" if empty? deref all then setref(all,reallyEnumerate()) deref(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} <>= )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} <>= )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 [members 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} <>= )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 => 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 => (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} <>= --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. @ <<*>>= <> -- Compile order for the differential equation solver: -- oderf.spad odealg.spad nlode.spad nlinsol.spad riccati.spad -- kovacic.spad lodof.spad odeef.spad <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}