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