\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra pfr.spad} \author{Robert Sutor, Barry Trager} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{domain PFR PartialFraction} <>= )abbrev domain PFR PartialFraction ++ Author: Robert S. Sutor ++ Date Created: 1986 ++ Change History: ++ 05/20/91 BMT Converted to the new library ++ Basic Operations: (Field), (Algebra), ++ coerce, compactFraction, firstDenom, firstNumer, ++ nthFractionalTerm, numberOfFractionalTerms, padicallyExpand, ++ padicFraction, partialFraction, wholePart ++ Related Constructors: ++ Also See: ContinuedFraction ++ AMS Classifications: ++ Keywords: partial fraction, factorization, euclidean domain ++ References: ++ Description: ++ The domain \spadtype{PartialFraction} implements partial fractions ++ over a euclidean domain \spad{R}. This requirement on the ++ argument domain allows us to normalize the fractions. Of ++ particular interest are the 2 forms for these fractions. The ++ ``compact'' form has only one fractional term per prime in the ++ denominator, while the ``p-adic'' form expands each numerator ++ p-adically via the prime p in the denominator. For computational ++ efficiency, the compact form is used, though the p-adic form may ++ be gotten by calling the function \spadfunFrom{padicFraction}{PartialFraction}. For a ++ general euclidean domain, it is not known how to factor the ++ denominator. Thus the function \spadfunFrom{partialFraction}{PartialFraction} takes as its ++ second argument an element of \spadtype{Factored(R)}. PartialFraction(R: EuclideanDomain): Cat == Capsule where FRR ==> Factored R SUPR ==> SparseUnivariatePolynomial R Cat == Join(Field, Algebra R,CoercibleTo Fraction R) with coerce: Fraction FRR -> % ++ coerce(f) takes a fraction with numerator and denominator in ++ factored form and creates a partial fraction. It is ++ necessary for the parts to be factored because it is not ++ known in general how to factor elements of \spad{R} and ++ this is needed to decompose into partial fractions. compactFraction: % -> % ++ compactFraction(p) normalizes the partial fraction \spad{p} ++ to the compact representation. In this form, the partial ++ fraction has only one fractional term per prime in the ++ denominator. firstDenom: % -> FRR ++ firstDenom(p) extracts the denominator of the first fractional ++ term. This returns 1 if there is no fractional part (use ++ \spadfunFrom{wholePart}{PartialFraction} to get the whole part). firstNumer: % -> R ++ firstNumer(p) extracts the numerator of the first fractional ++ term. This returns 0 if there is no fractional part (use ++ \spadfunFrom{wholePart}{PartialFraction} to get the whole part). nthFractionalTerm: (%,Integer) -> % ++ nthFractionalTerm(p,n) extracts the nth fractional term from ++ the partial fraction \spad{p}. This returns 0 if the index ++ \spad{n} is out of range. numberOfFractionalTerms: % -> Integer ++ numberOfFractionalTerms(p) computes the number of fractional ++ terms in \spad{p}. This returns 0 if there is no fractional ++ part. padicallyExpand: (R,R) -> SUPR ++ padicallyExpand(p,x) is a utility function that expands ++ the second argument \spad{x} ``p-adically'' in ++ the first. padicFraction: % -> % ++ padicFraction(q) expands the fraction p-adically in the primes ++ \spad{p} in the denominator of \spad{q}. For example, ++ \spad{padicFraction(3/(2**2)) = 1/2 + 1/(2**2)}. ++ Use \spadfunFrom{compactFraction}{PartialFraction} to return to compact form. partialFraction: (R, FRR) -> % ++ partialFraction(numer,denom) is the main function for ++ constructing partial fractions. The second argument is the ++ denominator and should be factored. wholePart: % -> R ++ wholePart(p) extracts the whole part of the partial fraction ++ \spad{p}. Capsule == add -- some constructor assignments and macros Ex ==> OutputForm fTerm ==> Record(num: R, den: FRR) -- den should have -- unit = 1 and only -- 1 factor LfTerm ==> List Record(num: R, den: FRR) QR ==> Record(quotient: R, remainder: R) Rep := Record(whole:R, fract: LfTerm) import %before?: (FRR,FRR) -> Boolean from Foreign Builtin -- private function signatures copypf: % -> % multiplyFracTerms: (fTerm, fTerm) -> % normalizeFracTerm: fTerm -> % partialFractionNormalized: (R, FRR) -> % -- declarations a,b: % n: Integer r: R -- private function definitions copypf(a: %): % == [a.whole,copy a.fract]$% LessThan(s: fTerm, t: fTerm): Boolean == -- have to wait until FR has < operation %before?(s.den,t.den) multiplyFracTerms(s : fTerm, t : fTerm) == nthFactor(s.den,1) = nthFactor(t.den,1) => normalizeFracTerm([s.num * t.num, s.den * t.den]$fTerm) : Rep i : Union(Record(coef1: R, coef2: R),"failed") coefs : Record(coef1: R, coef2: R) i := extendedEuclidean(expand t.den, expand s.den,s.num * t.num) i case "failed" => error "PartialFraction: not in ideal" coefs := (i :: Record(coef1: R, coef2: R)) c : % := copypf 0$% d : % if coefs.coef2 ~= 0$R then c := normalizeFracTerm ([coefs.coef2, t.den]$fTerm) if coefs.coef1 ~= 0$R then d := normalizeFracTerm ([coefs.coef1, s.den]$fTerm) c.whole := c.whole + d.whole not (null d.fract) => c.fract := append(d.fract,c.fract) c normalizeFracTerm(s : fTerm) == -- makes sure num is "less than" den, whole may be non-zero qr : QR := divide(s.num, (expand s.den)) qr.remainder = 0$R => [qr.quotient, nil()$LfTerm] -- now verify num and den are coprime f : R := nthFactor(s.den,1) nexpon : Integer := nthExponent(s.den,1) expon : Integer := 0 q : QR := divide(qr.remainder, f) while (q.remainder = 0$R) and (expon < nexpon) repeat expon := expon + 1 qr.remainder := q.quotient q := divide(qr.remainder,f) expon = 0 => [qr.quotient,[[qr.remainder, s.den]$fTerm]$LfTerm] expon = nexpon => (qr.quotient + qr.remainder) :: % [qr.quotient,[[qr.remainder, nilFactor(f,nexpon-expon)]$fTerm]$LfTerm] partialFractionNormalized(nm: R, dn : FRR) == -- assume unit dn = 1 nm = 0$R => 0$% dn = 1$FRR => nm :: % qr : QR := divide(nm, expand dn) c : % := [0$R,[[qr.remainder, nilFactor(nthFactor(dn,1), nthExponent(dn,1))]$fTerm]$LfTerm] d : % for i in 2..numberOfFactors(dn) repeat d := [0$R,[[1$R,nilFactor(nthFactor(dn,i), nthExponent(dn,i))]$fTerm]$LfTerm] c := c * d (qr.quotient :: %) + c -- public function definitions padicFraction(a : %) == b: % := compactFraction a null b.fract => b l : LfTerm := nil f : R e,d: Integer for s in b.fract repeat e := nthExponent(s.den,1) e = 1 => l := cons(s,l) f := nthFactor(s.den,1) d := degree(sp := padicallyExpand(f,s.num)) while (sp ~= 0$SUPR) repeat l := cons([leadingCoefficient sp,nilFactor(f,e-d)]$fTerm, l) d := degree(sp := reductum sp) [b.whole, sort(LessThan,l)]$% compactFraction(a : %) == -- only one power for each distinct denom will remain 2 > # a.fract => a af : LfTerm := reverse a.fract bf : LfTerm := nil bw : R := a.whole b : % s : fTerm := [(first af).num,(first af).den]$fTerm f : R := nthFactor(s.den,1) e : Integer := nthExponent(s.den,1) for t in rest af repeat f = nthFactor(t.den,1) => s.num := s.num + (t.num * (f **$R ((e - nthExponent(t.den,1)) : NonNegativeInteger))) b := normalizeFracTerm s bw := bw + b.whole if not (null b.fract) then bf := cons(first b.fract,bf) s := [t.num, t.den]$fTerm f := nthFactor(s.den,1) e := nthExponent(s.den,1) b := normalizeFracTerm s [bw + b.whole,append(b.fract,bf)]$% 0 == [0$R, nil()$LfTerm] 1 == [1$R, nil()$LfTerm] characteristic == characteristic$R coerce(r): % == [r, nil()$LfTerm] coerce(n): % == [(n :: R), nil()$LfTerm] coerce(a): Fraction R == q : Fraction R := (a.whole :: Fraction R) for s in a.fract repeat q := q + (s.num / (expand s.den)) q coerce(q: Fraction FRR): % == u : R := (recip unit denom q):: R r1 : R := u * expand numer q partialFractionNormalized(r1, u * denom q) a exquo b == b = 0$% => "failed" b = 1$% => a br : Fraction R := inv (b :: Fraction R) a * partialFraction(numer br,(denom br) :: FRR) recip a == (1$% exquo a) firstDenom a == -- denominator of 1st fractional term null a.fract => 1$FRR (first a.fract).den firstNumer a == -- numerator of 1st fractional term null a.fract => 0$R (first a.fract).num numberOfFractionalTerms a == # a.fract nthFractionalTerm(a,n) == l : LfTerm := a.fract (n < 1) or (n > # l) => 0$% [0$R,[l.n]$LfTerm]$% wholePart a == a.whole partialFraction(nm: R, dn : FRR) == nm = 0$R => 0$% -- move inv unit of den to numerator u : R := unit dn u := (recip u) :: R partialFractionNormalized(u * nm,u * dn) padicallyExpand(p : R, r : R) == -- expands r as a sum of powers of p, with coefficients -- r = HornerEval(padicallyExpand(p,r),p) qr : QR := divide(r, p) qr.quotient = 0$R => qr.remainder :: SUPR (qr.remainder :: SUPR) + monomial(1$R,1$NonNegativeInteger)$SUPR * padicallyExpand(p,qr.quotient) a = b == a.whole ~= b.whole => false -- must verify this (null a.fract) => null b.fract => a.whole = b.whole false null b.fract => false -- oh, no! following is temporary (a :: Fraction R) = (b :: Fraction R) - a == l: LfTerm := nil for s in reverse a.fract repeat l := cons([- s.num,s.den]$fTerm,l) [- a.whole,l] r * a == r = 0$R => 0$% r = 1$R => a b : % := (r * a.whole) :: % c : % for s in reverse a.fract repeat c := normalizeFracTerm [r * s.num, s.den]$fTerm b.whole := b.whole + c.whole not (null c.fract) => b.fract := append(c.fract, b.fract) b n * a == (n :: R) * a a + b == compactFraction [a.whole + b.whole, sort(LessThan,append(a.fract,copy b.fract))]$% a * b == null a.fract => a.whole * b null b.fract => b.whole * a af : % := [0$R, a.fract]$% -- a - a.whole c: % := (a.whole * b) + (b.whole * af) for s in a.fract repeat for t in b.fract repeat c := c + multiplyFracTerms(s,t) c coerce(a): Ex == null a.fract => a.whole :: Ex l : List Ex if a.whole = 0 then l := nil else l := [a.whole :: Ex] for s in a.fract repeat s.den = 1$FRR => l := cons(s.num :: Ex, l) l := cons(s.num :: Ex / s.den :: Ex, l) # l = 1 => first l reduce("+", reverse l) @ \section{package PFRPAC PartialFractionPackage} <>= )abbrev package PFRPAC PartialFractionPackage ++ Author: Barry M. Trager ++ Date Created: 1992 ++ BasicOperations: ++ Related Constructors: PartialFraction ++ Also See: ++ AMS Classifications: ++ Keywords: partial fraction, factorization, euclidean domain ++ References: ++ Description: ++ The package \spadtype{PartialFractionPackage} gives an easier ++ to use interfact the domain \spadtype{PartialFraction}. ++ The user gives a fraction of polynomials, and a variable and ++ the package converts it to the proper datatype for the ++ \spadtype{PartialFraction} domain. PartialFractionPackage(R): Cat == Capsule where -- R : UniqueFactorizationDomain -- not yet supported R : Join(EuclideanDomain, CharacteristicZero) FPR ==> Fraction Polynomial R INDE ==> IndexedExponents Symbol PR ==> Polynomial R SUP ==> SparseUnivariatePolynomial Cat == with partialFraction: (FPR, Symbol) -> Any ++ partialFraction(rf, var) returns the partial fraction decomposition ++ of the rational function rf with respect to the variable var. partialFraction: (PR, Factored PR, Symbol) -> Any ++ partialFraction(num, facdenom, var) returns the partial fraction ++ decomposition of the rational function whose numerator is num and ++ whose factored denominator is facdenom with respect to the variable var. Capsule == add partialFraction(rf, v) == df := factor(denom rf)$MultivariateFactorize(Symbol, INDE,R,PR) partialFraction(numer rf, df, v) makeSup(p:Polynomial R, v:Symbol) : SparseUnivariatePolynomial FPR == up := univariate(p,v) map(#1::FPR,up)$UnivariatePolynomialCategoryFunctions2(PR, SUP PR, FPR, SUP FPR) partialFraction(p, facq, v) == up := UnivariatePolynomial(v, Fraction Polynomial R) fup := Factored up ffact : List(Record(irr:up,pow:Integer)) ffact:=[[makeSup(u.factor,v) pretend up,u.exponent] for u in factors facq] fcont:=makeSup(unit facq,v) pretend up nflist:fup := fcont*(*/[primeFactor(ff.irr,ff.pow) for ff in ffact]) pfup:=partialFraction(makeSup(p,v) pretend up, nflist)$PartialFraction(up) coerce(pfup)$AnyFunctions1(PartialFraction up) @ \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. @ <<*>>= <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}