/usr/share/axiom-20170501/src/algebra/UPDECOMP.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 | )abbrev package UPDECOMP UnivariatePolynomialDecompositionPackage
++ Author: Frederic Lehobey
++ Date Created: 17 June 1996
++ Date Last Updated: 4 June 1997
++ References:
++ [1] Peter Henrici, Automatic Computations with Power Series,
++ Journal of the Association for Computing Machinery, Volume 3, No. 1,
++ January 1956, 10-15
++ [2] Dexter Kozen and Susan Landau, Polynomial Decomposition
++ Algorithms, Journal of Symbolic Computation (1989) 7, 445-456
-- Decomposition would be speeded up (O(n log n) instead of O(n^2)) by
-- implementing the algorithm described in [3] based on [4] and [5].
++ [3] Joachim von zur Gathen, Functional Decomposition Polynomials:
++ the Tame Case, Journal of Symbolic Computation (1990) 9, 281-299
++ [4] R. P. Brent and H. T. Kung, Fast Algorithms for Manipulating
++ Formal Power Series, Journal of the Association for Computing
++ Machinery, Vol. 25, No. 4, October 1978, 581-595
++ [5] R. P. Brent, Multiple-Precision Zero-Finding Methods and the
++ Complexity of Elementary Function Evaluation, Analytic
++ Computational Complexity, J. F. Traub, Ed., Academic Press,
++ New York 1975, 151-176
++ Description:
++ UnivariatePolynomialDecompositionPackage implements
++ functional decomposition of univariate polynomial with coefficients
++ in an \spad{IntegralDomain} of \spad{CharacteristicZero}.
UnivariatePolynomialDecompositionPackage(R,UP) : SIG == CODE where
R : Join(IntegralDomain,CharacteristicZero)
UP : UnivariatePolynomialCategory(R)
N ==> NonNegativeInteger
LR ==> Record(left: UP, right: UP)
QR ==> Record(quotient: UP, remainder: UP)
SIG ==> with
monicRightFactorIfCan : (UP,N) -> Union(UP,"failed")
++ monicRightFactorIfCan(f,d) returns a candidate to be the
++ monic right factor (h in f = g o h) of degree d of a
++ functional decomposition of the polynomial f or
++ \spad{"failed"} if no such candidate.
rightFactorIfCan : (UP,N,R) -> Union(UP,"failed")
++ rightFactorIfCan(f,d,c) returns a candidate to be the
++ right factor (h in f = g o h) of degree d with leading
++ coefficient c of a functional decomposition of the
++ polynomial f or \spad{"failed"} if no such candidate.
leftFactorIfCan : (UP,UP) -> Union(UP,"failed")
++ leftFactorIfCan(f,h) returns the left factor (g in f = g o h)
++ of the functional decomposition of the polynomial f with
++ given h or \spad{"failed"} if g does not exist.
monicDecomposeIfCan : UP -> Union(LR,"failed")
++ monicDecomposeIfCan(f) returns a functional decomposition
++ of the monic polynomial f of "failed" if it has not found any.
monicCompleteDecompose : UP -> List UP
++ monicCompleteDecompose(f) returns a list of factors of f for
++ the functional decomposition (\spad{[ f1, ..., fn ]} means
++ f = f1 o ... o fn).
CODE ==> add
rightFactorIfCan(p,dq,lcq) ==
dp := degree p
zero? lcq =>
error "rightFactorIfCan: leading coefficient may not be zero"
(zero? dp) or (zero? dq) => "failed"
nc := dp exquo dq
nc case "failed" => "failed"
n := nc::N
s := subtractIfCan(dq,1)::N
lcp := leadingCoefficient p
q: UP := monomial(lcq,dq)
k: N
for k in 1..s repeat
c: R := 0
i: N
for i in 0..subtractIfCan(k,1)::N repeat
c := c+(k::R-(n::R+1)*(i::R))*
coefficient(q,subtractIfCan(dq,i)::N)*
coefficient(p,subtractIfCan(dp+i,k)::N)
cquo := c exquo ((k*n)::R*lcp)
cquo case "failed" => return "failed"
q := q+monomial(cquo::R,subtractIfCan(dq,k)::N)
q
monicRightFactorIfCan(p,dq) == rightFactorIfCan(p,dq,1$R)
import UnivariatePolynomialDivisionPackage(R,UP)
leftFactorIfCan(f,h) ==
g: UP := 0
zero? degree h => "failed"
for i in 0.. while not zero? f repeat
qrf := divideIfCan(f,h)
qrf case "failed" => return "failed"
qr := qrf :: QR
r := qr.remainder
not ground? r => return "failed"
g := g+monomial(ground(r),i)
f := qr.quotient
g
monicDecomposeIfCan f ==
df := degree f
zero? df => "failed"
for dh in 2..subtractIfCan(df,1)::N | zero?(df rem dh) repeat
h := monicRightFactorIfCan(f,dh)
h case UP =>
g := leftFactorIfCan(f,h::UP)
g case UP => return [g::UP,h::UP]
"failed"
monicCompleteDecompose f ==
cf := monicDecomposeIfCan f
cf case "failed" => [ f ]
lr := cf :: LR
append(monicCompleteDecompose lr.left,[lr.right])
|