1)abbrev domain PFR PartialFraction 2++ Author: Robert S. Sutor 3++ Date Created: 1986 4++ Related Constructors: 5++ Also See: ContinuedFraction 6++ AMS Classifications: 7++ Keywords: partial fraction, factorization, euclidean domain 8++ References: 9++ Description: 10++ The domain \spadtype{PartialFraction} implements partial fractions 11++ over a euclidean domain \spad{R}. This requirement on the 12++ argument domain allows us to normalize the fractions. Of 13++ particular interest are the 2 forms for these fractions. The 14++ ``compact'' form has only one fractional term per prime in the 15++ denominator, while the ``p-adic'' form expands each numerator 16++ p-adically via the prime p in the denominator. For computational 17++ efficiency, the compact form is used, though the p-adic form may 18++ be gotten by calling the function 19++ \spadfunFrom{padicFraction}{PartialFraction}. For a general euclidean 20++ domain, it is not known how to factor the denominator. 21++ Thus the function \spadfunFrom{partialFraction}{PartialFraction} takes 22++ an element of \spadtype{Factored(R)} as its second argument. 23 24PartialFraction(R : EuclideanDomain) : Cat == Capsule where 25 FRR ==> Factored R 26 SUPR ==> SparseUnivariatePolynomial R 27 fTerm ==> Record(num : R, den : FRR) -- den should have unit = 1 28 -- and only 1 factor 29 LfTerm ==> List fTerm 30 NNI ==> NonNegativeInteger 31 32 Cat == Join(Field, Algebra R) with 33 coerce : % -> Fraction R 34 ++ coerce(p) sums up the components of the partial fraction and 35 ++ returns a single fraction. 36 37 coerce : Fraction FRR -> % 38 ++ coerce(f) takes a fraction with numerator and denominator in 39 ++ factored form and creates a partial fraction. It is 40 ++ necessary for the parts to be factored because it is not 41 ++ known in general how to factor elements of \spad{R} and 42 ++ this is needed to decompose into partial fractions. 43 44 compactFraction : % -> % 45 ++ compactFraction(p) normalizes the partial fraction \spad{p} 46 ++ to the compact representation. In this form, the partial 47 ++ fraction has only one fractional term per prime in the 48 ++ denominator. 49 50 numberOfFractionalTerms : % -> Integer 51 ++ numberOfFractionalTerms(p) computes the number of fractional 52 ++ terms in \spad{p}. This returns 0 if there is no fractional 53 ++ part. 54 55 fractionalTerms : % -> LfTerm 56 ++ fractionalTerms(p) extracts the fractional part of \spad{p} 57 ++ to a list of Record(num : R, den : Factored R). This returns 58 ++ [] if there is no fractional part. 59 60 padicallyExpand : (R, R) -> SUPR 61 ++ padicallyExpand(p, x) is a utility function that expands 62 ++ the second argument \spad{x} ``p-adically'' in 63 ++ the first. 64 65 padicFraction : % -> % 66 ++ padicFraction(q) expands the fraction p-adically in the primes 67 ++ \spad{p} in the denominator of \spad{q}. For example, 68 ++ \spad{padicFraction(3/(2^2)) = 1/2 + 1/(2^2)}. 69 ++ Use \spadfunFrom{compactFraction}{PartialFraction} to return to compact form. 70 71 partialFraction : (R, FRR) -> % 72 ++ partialFraction(numer, denom) is the main function for 73 ++ constructing partial fractions. The second argument is the 74 ++ denominator and should be factored. 75 76 wholePart : % -> R 77 ++ wholePart(p) extracts the whole part of the partial fraction 78 ++ \spad{p}. 79 80 if R has UniqueFactorizationDomain then 81 partialFraction : Fraction R -> % 82 ++ partialFraction(f) is a user friendly interface for partial 83 ++ fractions when f is a fraction of UniqueFactorizationDomain. 84 85 Capsule == add 86 87 -- some constructor assignments and macros 88 89 Ex ==> OutputForm 90 QR ==> Record(quotient : R, remainder : R) 91 92 Rep := Record(whole : R, fract : LfTerm) 93 94 -- private function signatures 95 96 copypf : % -> % 97 LessThan : (fTerm, fTerm) -> Boolean 98 multiplyFracTerms : (fTerm, fTerm) -> % 99 normalizeFracTerm : fTerm -> % 100 partialFractionNormalized : (R, FRR) -> % 101 102 -- private function definitions 103 104 -- define firstFactor/firstExponent because den of fTerm has only 1 factor 105 firstFactor(x : FRR) : R == first(factorList(x)).factor 106 firstExponent(x : FRR) : NNI == first(factorList(x)).exponent 107 108 copypf(a : %) : % == [a.whole, copy a.fract]$Rep 109 110 LessThan(s : fTerm, t : fTerm) == 111 -- have to wait until FR has < operation 112 if (GGREATERP(s.den, t.den)$Lisp pretend Boolean) then false 113 else true 114 115 multiplyFracTerms(s : fTerm, t : fTerm) == 116 firstFactor(s.den) = firstFactor(t.den) => 117 normalizeFracTerm([s.num * t.num, s.den * t.den]$fTerm) :: Rep 118 i : Union(Record(coef1 : R, coef2 : R),"failed") 119 coefs : Record(coef1 : R, coef2 : R) 120 i := extendedEuclidean(expand t.den, expand s.den, s.num * t.num) 121 i case "failed" => error "PartialFraction: not in ideal" 122 coefs := (i :: Record(coef1 : R, coef2 : R)) 123 c : % := copypf 0$% 124 d : % 125 if coefs.coef2 ~= 0$R then 126 c := normalizeFracTerm ([coefs.coef2, t.den]$fTerm) 127 if coefs.coef1 ~= 0$R then 128 d := normalizeFracTerm ([coefs.coef1, s.den]$fTerm) 129 c.whole := c.whole + d.whole 130 not(empty?(d.fract)) => c.fract := append(d.fract, c.fract) 131 c 132 133 normalizeFracTerm(s : fTerm) == 134 -- makes sure num is "less than" den, whole may be non-zero 135 qr : QR := divide(s.num, (expand s.den)) 136 qr.remainder = 0$R => [qr.quotient, []$LfTerm] 137 -- now verify num and den are coprime 138 f : R := firstFactor(s.den) 139 nexpon : NNI := firstExponent(s.den) 140 expon : NNI := 0 141 q : QR := divide(qr.remainder, f) 142 while (q.remainder = 0$R) and (expon < nexpon) repeat 143 expon := expon + 1 144 qr.remainder := q.quotient 145 q := divide(qr.remainder, f) 146 expon = 0 => [qr.quotient, [[qr.remainder, s.den]$fTerm]$LfTerm] 147 expon = nexpon => (qr.quotient + qr.remainder) :: % 148 [qr.quotient, [[qr.remainder, nilFactor(f, (nexpon-expon)::NNI)]$fTerm]$LfTerm] 149 150 partialFractionNormalized(nm : R, dn : FRR) == 151 -- assume unit dn = 1 152 nm = 0$R => 0$% 153 dn = 1$FRR => nm :: % 154 qr : QR := divide(nm, expand dn) 155 c : % := [0$R, [[qr.remainder, 156 nilFactor(firstFactor(dn), firstExponent(dn))]$fTerm]$LfTerm] 157 d : % 158 for i in rest factorList dn repeat 159 d := 160 [0$R, [[1$R, nilFactor(i.factor, i.exponent)]$fTerm]$LfTerm] 161 c := c * d 162 (qr.quotient :: %) + c 163 164 -- public function definitions 165 166 padicFraction(a : %) == 167 b : % := compactFraction a 168 empty?(b.fract) => b 169 l : LfTerm := [] 170 s : fTerm 171 f : R 172 e, d : Integer 173 for s in b.fract repeat 174 e := firstExponent(s.den) 175 e = 1 => l := cons(s, l) 176 f := firstFactor(s.den) 177 d := degree(sp := padicallyExpand(f, s.num)) 178 while (sp ~= 0$SUPR) repeat 179 l := cons([leadingCoefficient sp, nilFactor(f, (e-d)::NNI)]$fTerm, l) 180 d := degree(sp := reductum sp) 181 [b.whole, sort(LessThan, l)]$Rep 182 183 compactFraction(a : %) == 184 -- only one power for each distinct denom will remain 185 2 > # a.fract => a 186 af : LfTerm := reverse a.fract 187 bf : LfTerm := [] 188 bw : R := a.whole 189 b : % 190 s : fTerm := [(first af).num, (first af).den]$fTerm 191 f : R := firstFactor(s.den) 192 e : Integer := firstExponent(s.den) 193 t : fTerm 194 for t in rest af repeat 195 f = firstFactor(t.den) => 196 s.num := s.num + (t.num * 197 (f ^$R ((e - firstExponent(t.den)) :: NNI))) 198 b := normalizeFracTerm s 199 bw := bw + b.whole 200 if not(empty?(b.fract)) then bf := cons(first b.fract, bf) 201 s := [t.num, t.den]$fTerm 202 f := firstFactor(s.den) 203 e := firstExponent(s.den) 204 b := normalizeFracTerm s 205 [bw + b.whole, append(b.fract, bf)]$Rep 206 207 0 == [0$R, []$LfTerm] 208 1 == [1$R, []$LfTerm] 209 characteristic() == characteristic()$R 210 211 coerce(r : R) : % == [r, []$LfTerm] 212 coerce(n : Integer) : % == [(n :: R), []$LfTerm] 213 coerce(a : %) : Fraction R == 214 q : Fraction R := (a.whole :: Fraction R) 215 s : fTerm 216 for s in a.fract repeat 217 q := q + (s.num / (expand s.den)) 218 q 219 coerce(q : Fraction FRR) : % == 220 u : R := (recip unit denom q):: R 221 r1 : R := u * expand numer q 222 partialFractionNormalized(r1, u * denom q) 223 224 a : % exquo b : % == 225 b = 0$% => "failed" 226 b = 1$% => a 227 br : Fraction R := inv (b :: Fraction R) 228 a * partialFraction(numer br, (denom br) :: FRR) 229 recip a == (1$% exquo a) 230 231 numberOfFractionalTerms a == # a.fract 232 wholePart a == a.whole 233 234 fractionalTerms a == a.fract 235 236 partialFraction(nm : R, dn : FRR) == 237 nm = 0$R => 0$% 238 -- move inv unit of den to numerator 239 u : R := unit dn 240 u := (recip u) :: R 241 partialFractionNormalized(u * nm, u * dn) 242 243 padicallyExpand(p : R, r : R) == 244 -- expands r as a sum of powers of p, with coefficients 245 -- r = HornerEval(padicallyExpand(p, r), p) 246 qr : QR := divide(r, p) 247 qr.quotient = 0$R => qr.remainder :: SUPR 248 (qr.remainder :: SUPR) + monomial(1$R, 1$NonNegativeInteger)$SUPR * 249 padicallyExpand(p, qr.quotient) 250 251 a = b == 252 a.whole ~= b.whole => false -- must verify this 253 empty?(a.fract) => 254 empty?(b.fract) => a.whole = b.whole 255 false 256 empty?(b.fract) => false 257 -- oh, no! following is temporary 258 (a :: Fraction R) = (b :: Fraction R) 259 260 - a == 261 s : fTerm 262 l : LfTerm := [] 263 for s in reverse a.fract repeat l := cons([- s.num, s.den]$fTerm, l) 264 [- a.whole, l] 265 266 r : R * a : % == 267 r = 0$R => 0$% 268 r = 1$R => a 269 b : % := (r * a.whole) :: % 270 c : % 271 s : fTerm 272 for s in reverse a.fract repeat 273 c := normalizeFracTerm [r * s.num, s.den]$fTerm 274 b.whole := b.whole + c.whole 275 not(empty?(c.fract)) => b.fract := append(c.fract, b.fract) 276 b 277 278 n : Integer * a : % == (n :: R) * a 279 280 a + b == 281 compactFraction 282 [a.whole + b.whole, 283 sort(LessThan, append(a.fract, copy b.fract))]$Rep 284 285 a : % * b : % == 286 empty?(a.fract) => a.whole * b 287 empty?(b.fract) => b.whole * a 288 af : % := [0$R, a.fract]$Rep -- a - a.whole 289 c : % := (a.whole * b) + (b.whole * af) 290 s, t : fTerm 291 for s in a.fract repeat 292 for t in b.fract repeat 293 c := c + multiplyFracTerms(s, t) 294 c 295 296 coerce(a : %) : Ex == 297 empty?(a.fract) => a.whole :: Ex 298 s : fTerm 299 l : List Ex 300 if a.whole = 0 then l := [] else l := [a.whole :: Ex] 301 for s in a.fract repeat 302 s.den = 1$FRR => l := cons(s.num :: Ex, l) 303 l := cons(s.num :: Ex / s.den :: Ex, l) 304 # l = 1 => first l 305 reduce("+", reverse l) 306 307 if R has UniqueFactorizationDomain then 308 partialFraction f == partialFraction(numer f, factor denom f) 309 310)abbrev package PFRPAC PartialFractionPackage 311++ Author: Barry M. Trager 312++ Date Created: 1992 313++ BasicOperations: 314++ Related Constructors: PartialFraction 315++ Also See: 316++ AMS Classifications: 317++ Keywords: partial fraction, factorization, euclidean domain 318++ References: 319++ Description: 320++ The package \spadtype{PartialFractionPackage} gives an easier 321++ to use interface to \spadtype{PartialFraction}. 322++ The user gives a fraction of polynomials, and a variable and 323++ the package converts it to the proper datatype for the 324++ \spadtype{PartialFraction} domain. 325 326PartialFractionPackage(R) : Cat == Capsule where 327 R : Join(EuclideanDomain, PolynomialFactorizationExplicit, 328 CharacteristicZero) 329 INDE ==> IndexedExponents Symbol 330 PR ==> Polynomial R 331 FPR ==> Fraction PR 332 SUP ==> SparseUnivariatePolynomial 333 Cat == with 334 partialFraction : (FPR, Symbol) -> Any 335 ++ partialFraction(rf, var) returns the partial fraction decomposition 336 ++ of the rational function rf with respect to the variable var. 337 partialFraction : (PR, Factored PR, Symbol) -> Any 338 ++ partialFraction(num, facdenom, var) returns the partial fraction 339 ++ decomposition of the rational function whose numerator is num and 340 ++ whose factored denominator is facdenom with respect to the variable var. 341 Capsule == add 342 partialFraction(rf, v) == 343 df := factor(denom rf)$MultivariateFactorize(Symbol, INDE, R, PR) 344 partialFraction(numer rf, df, v) 345 346 makeSup(p : PR, v : Symbol) : SparseUnivariatePolynomial FPR == 347 up := univariate(p, v) 348 map((z1 : PR) : FPR +-> z1::FPR, up 349 )$UnivariatePolynomialCategoryFunctions2(PR, SUP PR, FPR, SUP FPR) 350 351 partialFraction(p, facq, v) == 352 up := UnivariatePolynomial(v, FPR) 353 fup := Factored up 354 fcont := makeSup(unit facq, v) pretend up 355 nflist : fup := fcont*(*/[primeFactor(makeSup(u.factor, v) pretend up,_ 356 u.exponent) for u in factorList facq]) 357 pfup := partialFraction(makeSup(p, v) pretend up, nflist 358 )$PartialFraction(up) 359 coerce(pfup)$AnyFunctions1(PartialFraction up) 360 361--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. 362--All rights reserved. 363-- 364--Redistribution and use in source and binary forms, with or without 365--modification, are permitted provided that the following conditions are 366--met: 367-- 368-- - Redistributions of source code must retain the above copyright 369-- notice, this list of conditions and the following disclaimer. 370-- 371-- - Redistributions in binary form must reproduce the above copyright 372-- notice, this list of conditions and the following disclaimer in 373-- the documentation and/or other materials provided with the 374-- distribution. 375-- 376-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the 377-- names of its contributors may be used to endorse or promote products 378-- derived from this software without specific prior written permission. 379-- 380--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 381--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 382--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 383--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 384--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 385--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 386--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 387--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 388--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 389--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 390--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 391