diff options
Diffstat (limited to 'src/algebra/fparfrac.spad.pamphlet')
| -rw-r--r-- | src/algebra/fparfrac.spad.pamphlet | 232 |
1 files changed, 232 insertions, 0 deletions
diff --git a/src/algebra/fparfrac.spad.pamphlet b/src/algebra/fparfrac.spad.pamphlet new file mode 100644 index 00000000..9afe1d78 --- /dev/null +++ b/src/algebra/fparfrac.spad.pamphlet @@ -0,0 +1,232 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra fparfrac.spad} +\author{Manuel Bronstein} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{domain FPARFRAC FullPartialFractionExpansion} +<<domain FPARFRAC FullPartialFractionExpansion>>= +)abbrev domain FPARFRAC FullPartialFractionExpansion +++ Full partial fraction expansion of rational functions +++ Author: Manuel Bronstein +++ Date Created: 9 December 1992 +++ Date Last Updated: 6 October 1993 +++ References: M.Bronstein & B.Salvy, +++ Full Partial Fraction Decomposition of Rational Functions, +++ in Proceedings of ISSAC'93, Kiev, ACM Press. +FullPartialFractionExpansion(F, UP): Exports == Implementation where + F : Join(Field, CharacteristicZero) + UP : UnivariatePolynomialCategory F + + N ==> NonNegativeInteger + Q ==> Fraction Integer + O ==> OutputForm + RF ==> Fraction UP + SUP ==> SparseUnivariatePolynomial RF + REC ==> Record(exponent: N, center: UP, num: UP) + ODV ==> OrderlyDifferentialVariable Symbol + ODP ==> OrderlyDifferentialPolynomial UP + ODF ==> Fraction ODP + FPF ==> Record(polyPart: UP, fracPart: List REC) + + Exports ==> Join(SetCategory, ConvertibleTo RF) with + "+": (UP, $) -> $ + ++ p + x returns the sum of p and x + fullPartialFraction: RF -> $ + ++ fullPartialFraction(f) returns \spad{[p, [[j, Dj, Hj]...]]} such that + ++ \spad{f = p(x) + \sum_{[j,Dj,Hj] in l} \sum_{Dj(a)=0} Hj(a)/(x - a)\^j}. + polyPart: $ -> UP + ++ polyPart(f) returns the polynomial part of f. + fracPart: $ -> List REC + ++ fracPart(f) returns the list of summands of the fractional part of f. + construct: List REC -> $ + ++ construct(l) is the inverse of fracPart. + differentiate: $ -> $ + ++ differentiate(f) returns the derivative of f. + D: $ -> $ + ++ D(f) returns the derivative of f. + differentiate: ($, N) -> $ + ++ differentiate(f, n) returns the n-th derivative of f. + D: ($, NonNegativeInteger) -> $ + ++ D(f, n) returns the n-th derivative of f. + + Implementation ==> add + Rep := FPF + + fullParFrac: (UP, UP, UP, N) -> List REC + outputexp : (O, N) -> O + output : (N, UP, UP) -> O + REC2RF : (UP, UP, N) -> RF + UP2SUP : UP -> SUP + diffrec : REC -> REC + FP2O : List REC -> O + +-- create a differential variable + u := new()$Symbol + u0 := makeVariable(u, 0)$ODV + alpha := u::O + x := monomial(1, 1)$UP + xx := x::O + zr := (0$N)::O + + construct l == [0, l] + D r == differentiate r + D(r, n) == differentiate(r,n) + polyPart f == f.polyPart + fracPart f == f.fracPart + p:UP + f:$ == [p + polyPart f, fracPart f] + + differentiate f == + differentiate(polyPart f) + construct [diffrec rec for rec in fracPart f] + + differentiate(r, n) == + for i in 1..n repeat r := differentiate r + r + +-- diffrec(sum_{rec.center(a) = 0} rec.num(a) / (x - a)^e) = +-- sum_{rec.center(a) = 0} -e rec.num(a) / (x - a)^{e+1} +-- where e = rec.exponent + diffrec rec == + e := rec.exponent + [e + 1, rec.center, - e * rec.num] + + convert(f:$):RF == + ans := polyPart(f)::RF + for rec in fracPart f repeat + ans := ans + REC2RF(rec.center, rec.num, rec.exponent) + ans + + UP2SUP p == + map(#1::UP::RF, p)$UnivariatePolynomialCategoryFunctions2(F, UP, RF, SUP) + + -- returns Trace_k^k(a) (h(a) / (x - a)^n) where d(a) = 0 + REC2RF(d, h, n) == +-- one?(m := degree d) => + ((m := degree d) = 1) => + a := - (leadingCoefficient reductum d) / (leadingCoefficient d) + h(a)::UP / (x - a::UP)**n + dd := UP2SUP d + hh := UP2SUP h + aa := monomial(1, 1)$SUP + p := (x::RF::SUP - aa)**n rem dd + rec := extendedEuclidean(p, dd, hh)::Record(coef1:SUP, coef2:SUP) + t := rec.coef1 -- we want Trace_k^k(a)(t) now + ans := coefficient(t, 0) + for i in 1..degree(d)-1 repeat + t := (t * aa) rem dd + ans := ans + coefficient(t, i) + ans + + fullPartialFraction f == + qr := divide(numer f, d := denom f) + qr.quotient + construct concat + [fullParFrac(qr.remainder, d, rec.factor, rec.exponent::N) + for rec in factors squareFree denom f] + + fullParFrac(a, d, q, n) == + ans:List REC := empty() + em := e := d quo (q ** n) + rec := extendedEuclidean(e, q, 1)::Record(coef1:UP,coef2:UP) + bm := b := rec.coef1 -- b = inverse of e modulo q + lvar:List(ODV) := [u0] + um := 1::ODP + un := (u1 := u0::ODP)**n + lval:List(UP) := [q1 := q := differentiate(q0 := q)] + h:ODF := a::ODP / (e * un) + rec := extendedEuclidean(q1, q0, 1)::Record(coef1:UP,coef2:UP) + c := rec.coef1 -- c = inverse of q' modulo q + cm := 1::UP + cn := (c ** n) rem q0 + for m in 1..n repeat + p := retract(em * un * um * h)@ODP + pp := retract(eval(p, lvar, lval))@UP + h := inv(m::Q) * differentiate h + q := differentiate q + lvar := concat(makeVariable(u, m), lvar) + lval := concat(inv((m+1)::F) * q, lval) + qq := q0 quo gcd(pp, q0) -- new center + if (degree(qq) > 0) then + ans := concat([(n + 1 - m)::N, qq, (pp * bm * cn * cm) rem qq], ans) + cm := (c * cm) rem q0 -- cm = c**m modulo q now + um := u1 * um -- um = u**m now + em := e * em -- em = e**{m+1} now + bm := (b * bm) rem q0 -- bm = b**{m+1} modulo q now + ans + + coerce(f:$):O == + ans := FP2O(l := fracPart f) + zero?(p := polyPart f) => + empty? l => (0$N)::O + ans + p::O + ans + + FP2O l == + empty? l => empty() + rec := first l + ans := output(rec.exponent, rec.center, rec.num) + for rec in rest l repeat + ans := ans + output(rec.exponent, rec.center, rec.num) + ans + + output(n, d, h) == +-- one? degree d => + (degree d) = 1 => + a := - leadingCoefficient(reductum d) / leadingCoefficient(d) + h(a)::O / outputexp((x - a::UP)::O, n) + sum(outputForm(makeSUP h, alpha) / outputexp(xx - alpha, n), + outputForm(makeSUP d, alpha) = zr) + + outputexp(f, n) == +-- one? n => f + (n = 1) => f + f ** (n::O) + +@ +\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>> + +<<domain FPARFRAC FullPartialFractionExpansion>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |