FAST ALGORITHMS FOR PARTIAL FRACTION DECOMPOSITION

The partial fraction decomposition of a proper rational function whose denominator has degree n and is given in general factored form can be done in 0(n log n) operations in the worst case. 2 require 0(n ) operations, and 0(n log n) operations for the special case where the factors appearing in the denominator are all linear.

The general partial fraction compute the coefficients of the polynomials C for i j ° 1 9 • • •, SL^ such that I-em) is to 1,...,k and tr-l Z is i n Q?(x) i-i j-i Q i (x) with deg C < deg Q for all i,j. The existence and uniqueness of the polypi J nomials C. . ate well known (see, e.g., van der Waerden [1953]). There are i* J enormous applications of partial fractions in applied mathematics and in network theory (see, e.g., Henrici [1974] and Weinberg [1962]). This paper gives fast algorithms for solving the general partial fraction expansion problem when n is large.
Recently Chin and Ullman [1975] showed that in case that all the Q i 3/2 have degree one the problem can be done in 0((n log n) ) operations. This bound was further improved by Chin in his thesis (Chin [1975]). He showed that if the Q i are all linear, then the problem can be done in O((log k)-(n log n)) operations. However, the assumption that the are all linear factors is crucial in his methods. Hence the problem of solving the general PF problem 2 (without assuming that the Q i are linear) in 0(n ) operations is stated as an unsolved problem in his thesis. Note that the general PF problem does occur frequently in practice. For example, if we work over the field of real numbers, then the factors Q i certainly can have either degree one or two. (See also Grau [1971] and Henrici [1971] for more examples.) In this paper, we show that the general PF problem can done in O((log n)»M(n)) operations in the worst case. M(n) is any upper bound on the number of operations needed to multiply two nth degree polynomials, which satisfies some mild regularity condition (see Section 2). In particular, if an FFT algorithm is used for polynomial multiplication (see, e.g., Knuth [1969], Borodin and Munro [1975]), then we have M(n) = 0(n log n), which satisfies the regularity condition, and 2 hence the general PF problem can be done in 0(n log n) operations. Moreover, we note that for the special case where the are all linear, our approach will lead to Chin's O((log k)•(n log n)) algorithm.
Basic assumptions and preliminary lemmas used in this paper are introduced in Section 2. In Section 3, the solution of the general PF problem is reduced to the solution of two simpler problems, problem PI and problem P2, - Let Z* be the set of all nonnegative integers and let G: Z + -* Z + be a nondecreasing function. We say G satisfies Condition C» if G(n) -n.H(n) for some nondecreasing function H: Z + Z + . We assume that M satisfies Condition C. Similar regularity conditions are usually assumed (see, e.g., Aho, Hopcroft and Ullman [1974, p.280], Brent and Kung [1976] and Moenck [1973b]).
There are many algorithms for polynomial multiplication. For example, the 2 classical algorithm gives M(n) » c^n , binary splitting multiplication gives 1 585 M(n) « c 2 n , and if the field K is algebraically closed, then FFT mult* plication gives M(n) « c 3 n log n, where c.j, c 2 , c 3 are positive constants (see e.g., Fateman [1974]). In all cases M satisfies Condition C. In fact all we need in this paper are some consequences of Condition C. Hence it is possible to weaken our assumption on M, if one wishes to do so.
-5-Let D(n) be the number of operations needed to divide a polynomial of degree 2n by a polynomial of degree n. Then using Newton's method and the fact that M satisfies Condition C, one can show the following lemma (see, e.g., Borodin and Munro [1975] and Kung [1974]): Using the algorithm EGCD in Moenck [1973a], which is a generalization of an algorithm due to Schoenhage [1971] for integer GCDs, one can show the following lemma.
We shall assume that F satisfies the condition that for any n ± £ Z + . Clearly, if F(n) « c-(log n)*M(n) for some positive constant c as in Lemma 2.2, then F satisfies the condition. In fact, the required condition in F is satisfied as long as F satisfies Condition C.

PROBLEMS Pi, P2 AND STATEMENT OF RESULTS
Consider the following two instances of the general PF problem defined in Section 1.
Problem PI : (This is the general PF problem with JL » 1 for all i.) k Given the fraction P/ Tl R. where the R. are relatively prime and The decomposition (3.1) is called the incomplete partial fraction decomposition by Henrici [1971Henrici [ , 1974. Note also that efficient algorithms for solving problem PI will furnish efficient procedures for factoring polynomials, as observed by Grau [1971].
Problem P2? (This is the general PF problem with k = 1.) Given the fraction P/Q where deg P < jfc»deg Q compute the polynomials C|,... ,C^ such that The following lemma essentially shows that fast algorithms for problems PI and P2 will lead to fast algorithms for the general PF problem. Define T(k,n),T 1 (k,n) and T (£,deg Q) to be the number of operations needed to solve The result follows from the observation that general PF problem can be solved in the following way: 1. Multiply Q. (x) out for i -1,...,k. Let the expansion of (^(x) be R j ,(x) for all i.
Note that each Q ± (x) can be computed in (KM^-deg Q i >) operations by an algorithm in Brent [1975].
• We summarize our results on T^k.n) and T 2 U,deg Q) in the following: We have the following results for the general partial fraction expansion problem.
Theorem 3.1 The general PF problem can be done in The result follows from (i), (iii) and Lemma 3.1. • Corollary 3 J 2 The general PF problem can be done in 0(n log n) operations.

Proof
Note that in Theorem 3.1, k ^ n and & £ n. The result follows from the theorem and Lemma 2.2 by letting M(n) • 0(n log n).
• 2 0(n log n) is the best asymptotic bound known for the general PF problem.
Theorem 3.2 The general PF problem can be done in O((log k)•(n log n)) operations, if Q i (x) -x-z i , for i»l,...,k. -9-

Proof
The result follows from (ii), ( iv ) and Lemma 3.1.
The bound in Theorem 3.2 was obtained previously by Chin [1975]. We include it here just to show that his result will emerge as a special case in our general approach. See the remarks at the end of Section 4. -10-

AN ALGORITHM FOR PROBLEM Pi
We first assume that P(x) « 1 in problem Pi. Thus we want to find   On the other hand, suppose that deg ^ 1. Then by (4.6), (4.7) and (4.8), we know that the coefficients of A, and those of satisfy the same system of equations. Because the R^ are relatively prime, the system is nonsingular. In the following we study the number of operations needed by the algorithm. k It is well known that II R.(x) can be computed by using a binary split-i=1 1 ting scheme, which is illustrated as follows for the case k » 8:

Proof
We shall again use the binary splitting technique. We may assume that k is a power of 2. It is easy to check that This gives us a recursive algorithm for computing R. By Lemma 4.1, we may assume that all the products such as II R and n R needed by the j>k/2 J J u algorithm have been precomputed. The result again follows from the fact that the of the degrees of all polynomials at any level of the associated binary sum tree is n. Hence we have the following algorithm for computing A^: Algorithm 4.2 1. Compute D i (x) such that D^x) « D^x+z^.
A 3. Compute A^(x) such that A^x) -A^x-z^.
Step 1 is equivalent to evaluating D^^ and all its derivatives at z^. Aho, Steiglitz and Ullman [1975] and Vari [1974] have independently shown that this can be done in 0(m^ log m^) operations. Similarly, step 3 can be done in 0(m. log m.) operations.
Step 2 involves a division, which is 0(m log m )   operations.

A SPECIAL CASE FOR PROBLEM P2
The following theorem can be found in Chin and Ullman [1975], Theorem 6.1 Problem P2 can be solved In 0(jj log 1) operations if deg Q « 1, In this section we extend the theorem to the case that deg Q « 2, Our result is of interest when the underlying field K is the field of real numbers, for in this case irreducible factors can have either degree one or two. We may assume that Q is monic, since this will effect only 0(A) operations.
Step 3 clearly uses 0(1) operations. Thus, we have shown the following theorem.
Theorem 6.2 Problem P2 can be solved in 0(11 log Z) operations if deg Q « 2.
It is an open problem whether Theorem 6.2 holds if deg Q > 2.