/usr/share/axiom-20170501/src/algebra/LODOF.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 | )abbrev package LODOF LinearOrdinaryDifferentialOperatorFactorizer
++ Author: Fritz Schwarz, Manuel Bronstein
++ Date Created: 1988
++ Date Last Updated: 3 February 1994
++ References:
++ Bron95 On radical solutions of linear ordinary differential equations
++ Abra01 On Solutions of Linear Functional Systems
++ Muld95 Primitives: Orepoly and Lodo
++ Description:
++ \spadtype{LinearOrdinaryDifferentialOperatorFactorizer} provides a
++ factorizer for linear ordinary differential operators whose coefficients
++ are rational functions.
LinearOrdinaryDifferentialOperatorFactorizer(F, UP) : SIG == CODE where
F : Join(Field, CharacteristicZero,
RetractableTo Integer, RetractableTo Fraction Integer)
UP: UnivariatePolynomialCategory F
RF ==> Fraction UP
L ==> LinearOrdinaryDifferentialOperator1 RF
SIG ==> 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.
CODE ==> add
import RationalLODE(F, UP)
import RationalRicDE(F, UP)
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 ((leadingCoefficient r) = 1) =>
recurfactor(op, r, zeros, ezfactor, adj?)
op1 := opeval(op, dd - coefficient(r, 0)::L)
map_!((z:L):L+->opeval(z,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) ==
(n = 1) =>
(u := expsol(l, zeros, ezfactor)) case "failed" => "failed"
D() - u::RF::L
"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((z1:F):F+->z1,p)_
$UnivariatePolynomialCategoryFunctions2(F, UP,
F, SparseUnivariatePolynomial F))]
if F is AlgebraicNumber then
import AlgFactor UP
factor l ==
innerFactor(l,(p:UP):List(F)+->zro(p,factor),factor,true)
factor1 l ==
innerFactor(l,(p:UP):List(F)+->zro(p,factor),factor,false)
else
factor l ==
innerFactor(l,(p:UP):List(F)+->zro(p,squareFree),squareFree,true)
factor1 l ==
innerFactor(l,(p:UP):List(F)+->zro(p,squareFree),squareFree,false)
|