Proceedings of the Steklov Institute of Mathematics, Vol. 218, 1997, pp. 117–128. Translated from Trudy Matematicheskogo Instituta imeni V. A. Steklova, Vol. 218, 1997, pp. 122–133. Original Russian Text Copyright \begin{matrix}c \\ ◯ \end{matrix} 1997 by Gashkov. English Translation Copyright \begin{matrix}c \\ ◯ \end{matrix} by ���� �����/Interperiodica Publishing (Russia).

On the Complexity of Integration of Rational Fractions

S. B. Gashkov

Received February, 1997

1. Introduction. The algortithm of integration of rationl fractions is one of the most significant and most often used algorithms of analysis. It is well known that the integral of a rational fraction can be represented as the sum of a rational and a logarithmic part, and that the rational part can be separated only by means of arithmetic operations. The algorithm for this purpose was found by Ostrogradskii in 1844 (see [1]). Twenty years later Hermite proposed a different algorithm for the same purpose (see [3]). In the Russian literature the algorithm for separating the rational part of the integral of a rational fraction is called the Hermite–Ostrogradskii algorithm, but actually the Hermite algorithm is usually presented. In foreign literature the name of Ostrogradskii is usually omitted, although his article [1] was published in St. Petersbourg in French.

After Karatsuba proposed in 1962 a fast algorithm for multiplication of numbers [4] (in his recent article [5] he himself writes, in an interesting form, how it happened) an interest arose to the evaluation of the complexity of algorithms of algebra and analysis and to the consruction of fast modifications of these algorithms. Different fast algorithms have been found for the main arithmetic operations with numbers and polynomials in the works [6–10] and others, as well as fast algorithms for computing the values of polynomials and their interpolations (see, e.g., [11]), fast algorithms for computing determinants and for solving systems of linear equations (after [12] dozens of works were written on this subject), and many other results, among which there are also those which can be considered as polynomial algorithms. It is impossible to cite all of them in a short article (see, e.g., [13–15] and the references therein).

It is shown in [16] how Hermite’s algorithm can be modernized so that it could operate with the complexity O(n log^{2} n), where n is the degree of the denominator of the fraction being integrated (which can be assumed to be proper). This algorithm consists of four subalgorithms, namely, a square-free fractorization of the denominator of the fraction, its fractorization into partial fractions, integration of the obtained partial fractions by parts and the summation of the rational parts of the resulting integrals and integrands which are also rational fractions. Each of these algorithms is of interest in itself and has numerous applications, especially the first two and the last one. All these algorithms have the same estimate of the complexity in the order of magnitude.

In the work [16] only the estimates of the complexity of the first and the third algorithm are given in detail since the estimates for the other two algorithms were obtained in [17, 18]. These three works contain a long history of this question, and therefore we omit it here.

However, it is not written there that this history dates back to Ostrogradskii since the algorithm for square-free factorization [16] is given in [1] (the author of [16] evidently does not know this fact). This algorithm is somewhat more complicated for understanding than the standard algorithm (which dates back to the Dutch mathematician of the century before last Johann Hudde), but is, nevertheless, simple enough for its presentation to first-year students, and therefore deserves a wider propagation and application.

2. Estimate of the complexity of the algorithm for rational fractions integration and its subalgorithms. We can propose the estimate of the complexity of integration of regular rational fractions of degree n in the form

⌈log_{2} k⌉(21 M(n) + 11n) + 24 M(n) + 26n + 4 GD(n),

where k is the number of square-free factors in Ostrogradskii’s factorization, M(n) is the estimate of the complexity of multiplication of two polynomials of degree n satisfying the condition of superadditivity: M(x) + M(y) ≤ M(x + y), and GD(n) is a similar estiamte of complexity of computing the GCD of two polynomials. It is known (see, e.g., [13]) that we can take Cn log n as M(n) and C M(n) log n as GD(n) for a large C (this was proved in [19, 20]).

In this case, we can estimate the complexity of square-free factorization of the polynomial P of degree n as 2 GD(n) + 7 M(n) + 6n, the complexity of partial fraction decomposition of a regular fraction with denominator P as

GD(n) + ⌈log_{2} k⌉(17 M(n) + \frac{17n}{2}) + 8 M(n) + 5n,

the complxity of integration of the resulting partial fractions as

3 M(n) + 6 M(p) + GD(n − p) + 15p,

where p = deg(P, dP/dx), and the complexity of summation of the partial fractions resulting from this intergration as

⌈log_{k} t⌉(3 M(n) + M(p) + \frac{5p}{2}) + M(n − p) + 3(n − p).

3. Ostrogradskii’s algorithm for square-free factorization of polynomials. We shall give the algorithm for square-free factorization in the notation of Ostrogradskii and estimate its complexity (i.e., the number of operations it performs). Let P = P_{1} be a polynomial being factorized. We set

\begin{matrix}P_{2} = (P_{1}, \frac{dP_{1}}{dx} ),  Q_{1} = \frac{P_{1}}{P_{2}} ,  R_{1} = \frac{dP_{1}/dx}{P_{2}},  q_{1} = (Q_{1}, R_{1} − \frac{dQ_{1}}{dx} ), \\ Q_{2} = \frac{Q_{1}}{q_{1}} ,  R_{2} = \frac{R_{1} − dQ_{1}/dx}{q_{1}},  q_{2} = (Q_{2}, R_{2} − \frac{dQ_{2}}{dx} ) \end{matrix}

and so on, until we get R_{k+1} = dQ_{k+1}/dx (here and in the sequel we denote by (f, g) the GCD of the polynomials f and g). Then P = q_{1}q_{2}^{2} ⋯ q_{k}^{k} is the required decomposition since q_{i} do not have multiple roots and are coprime. Note that we have incidentally found the decomposition P = Q_{1}Q_{2} ⋯ Q_{k}, where Q_{i} = q_{i} ⋯ q_{k}, i = 1, …, k, are also square-free factors. However, they are no longer coprime but successively divide one another.

Let us compare this algorithm with the standard one. In the latter, the sequence of polynomials is recurrently computed from the relation P_{i+1} = (P_{i}, dP_{i}/dx) and, with the use of the theorem on the lowering of the multiplicity of the root in the process of differentiation, it is verified that

P_{i} = q_{i}q_{i+1}^{2} ⋯ q_{k}^{k−i+1} = Q_{i}Q_{i+1} ⋯ Q_{k},

and then it is found by means of division that

Q_{i} = \frac{P_{i}}{P_{i+1}} , q_{i} = \frac{Q_{i}}{Q_{i+1}} .

Ostrogradskii’s algorithm is substantiated by his remarkable identities

R_{i+1} = \frac{dP_{i+1}/dx}{P_{i+2}} = \frac{R_{i} − dQ_{i}/dx}{q_{i}}, q_{i} = (Q_{i}, R_{i} − \frac{dQ_{i}}{dx} ).

For the convenience of the reader we shall give their proof, the more so as it is short. Using induction, we find that

R_{i} − \frac{dQ_{i}}{dx} = \frac{dP_{i}/dx}{P_{i+1}} − \frac{d(P_{i}/P_{i+1})}{dx} = \frac{P_{i} × (dP_{i+1}/dx)}{P_{i+1}^{2}} = \frac{Q_{i} × (dP_{i+1}/dx)}{P_{i+1}} = \frac{Q_{i+1}q_{i} × (dP_{i+1}/dx)}{P_{i+1}},

whence we have

R_{i+1} = \frac{R_{i} − dQ_{i}/dx}{q_{i}} = \frac{Q_{i+1} × (dP_{i+1}/dx)}{P_{i+1}} = \frac{dP_{i+1}/dx}{P_{i+2}}.

Since

P_{i} = q_{i}q_{i+1}^{2} ⋯ q_{k}^{k−i+1},  R_{i} = \frac{dP_{i}/dx}{P_{i+1}},  Q_{i} = q_{i} ⋯ q_{k},

it follows from the theorem on the lowering of the multiplicity of the root in the process of differentiation that (R_{i}, Q_{i}) = 1, whence it follows, now without any induction, that

(Q_{i}, R_{i} − \frac{dQ_{i}}{dx} ) = (Q_{i}, R_{i+1}q_{i}) = (Q_{i+1}q_{i}, R_{i+1}q_{i}) = (Q_{i+1}q_{i}, R_{i+1}q_{i}) = q_{i}(Q_{i+1}, R_{i+1}) = q_{i}.

Note that the relations

R_{i} = \frac{dP_{i}/dx}{P_{i+1}}, Q_{i} = \frac{P_{i}}{P_{i+1}}

imply that deg R_{i} = deg P_{i} − deg P_{i+1} − 1 = deg Q_{i} − 1.

The advantage of Ostrogradskii’s algorithm as compared to the astandard one is obvious since at every step of iteration it requires the finding of the GCD of the lower-degree polynomials than in the standard algorithm, and one division is used rather than two, and again for polynomials of a lower degree. Morever, it admits of a good evaluation of the complexity when fast algorithms of multiplication, division, and calculation of the GCD are used in it.

4. Estimate of the complexity of Ostrogradskii’s algorithm. It can be obtined very simply. We set deg Q_{i} = n_{i}, deg q_{i} = m_{i}, deg P_{2} = p = n_{2} + ⋯ + n_{k}, and then m_{i} = n_{i} − n_{i+1}, deg R_{i} = n_{i} − 1. We estimate the complexity of the exact division of a polynomial of degree m by a polynomial of degree n as 3(M(m − n) + m − n) + M(m) (see, e.g., [21, problems 17.25, 17.32]), and, hence, the complexity of the division of two different polynomials of a degree not higher than m by the same polynomial of degree n can be estimated as 3(M(m − n) + m − n) + 2M(m) (see, e.g., [21, solutions of problems 17.25–17.26]). Then the required estimate has the form

n + GD(n) + 3 M(n − p) + 3n − 3p + 2 M(n) + ∑_{i=1}^{k}(GD(n_{i}) + 2 M(n_{i}) + 2n_{i}) + ∑_{i=2}^{k}(3 M(n_{i}) + 3n_{i})

≤ n + GD(n) + 3 M(n − p) + 3n − 3p + 2 M(n) + GD(n) + 2 M(n) + 2n + 3 M(p) + 3p

≤ 2 GD(n) + 7 M(n) + 6n.

The obtained estimate shows that the complexity of Ostrogradskii’s algorithm for a squarefree factorization of an nth degree polynomial exceeds only twice the complexity of the extended algorithm for finding the GCD of two polynomials of degree n.

5. One more version of Ostrogradskii’s algorithm. It is interesting to note that Ostogradskii propossed in his note [2] one more version of the factorization algorithm based on the identity q_{j} = (Q, R − j × (dQ/dx)), where

Q = Q_{1} = \frac{P_{1}}{P_{2}} = \frac{P}{P_{2}} ,  R = R_{1} = \frac{dP_{1}/dx}{P_{2}},  P_{2} = (\frac{dP_{1}}{dx} , P_{1}).

The proof follows from the chain of relations (in which, for brevity, the differentiation is denote by prime)

\begin{matrix}R = \frac{(P_{1})′}{P_{2}} = \frac{(q_{1}q_{2}^{2} ⋯ q_{k}^{k})′}{q_{2} ⋯ q_{k−1}^{k}} = ∑_{i=1}^{k} iq_{1} ⋯ q_{i−1}(q_{i})′q_{i+1} ⋯ q_{k}, \\ Q′ = ∑_{i=1}^{k} q_{1} ⋯ q_{i−1}(q_{i})′q_{i+1} ⋯ q_{k}, R − jQ′ = ∑_{i=1}^{k}(i − j)q_{1} ⋯ q_{i−1}(q_{i})′q_{i+1} ⋯ q_{k} \end{matrix}

(justified by the Leibniz differentiation rule) and from the note stating that all terms in the last sum, except for the sth term are divisible by q_{s}, and the sth term is prime to q_{s} (since all its factors are prime to it). Consequently, (R − jQ′, q_{s}) = 1 for s ≠ j and, obviously, the polynomial R − jQ is divisible by the polynomial q_{j}.

The indicated Ostogradskii formulas are more elegant than the previous formulas, but the algorithm described in this section is obviously more complicated than the preceding one since the degrees of the polynomials in the sequence in it, for which we have to find the GCD, do not become lower, and it does not admit of an obvious good upper bound. However, if we have to calculate not all factors q_{i}, but only any one of them, then this algorithm is obviously the fastest.

6. Estimate of the complexity of the algorithm for factorization of a rational fraction into partial fractions. The main idea of the fast algorithm for factorization into partial fractions consists in the application of Karatsuba’s method of division in half.

Let P/Q be an arbitrary proper fraction of degree n, Q = q_{1} ⋯ q_{s}, (q_{i}, q_{j}) = 1, i ≠ j. We set Q_{i} = Q/q_{i}, and then, as is known, the relation

\frac{P}{Q} = ∑_{i=1}^{s} \frac{p_{i}}{q_{i}} ,

where

p_{i} = (P mod q_{i} × Q_{i}^{−1} mod q_{i}) mod q_{i}.

holds true. In order to substantiate this relation, we must verify that (Q_{1}, …, Q_{s}) = 1, use the linear representation of the GCD 1 = ∑_{i=1}^{s} U_{i}Q_{i}, and note that U_{i} mod q_{i} = Q_{i}^{−1} mod q_{i} and

\frac{P}{Q} = ∑_{i=1}^{s} \frac{P × U_{i}}{q_{i}} = ∑_{i=1}^{s} \frac{(P × U_{i}) mod q_{i}}{q_{i}} .

Assume that we have already calculated P mod q_{i} and Q_{i} mod q_{i}. Then Q_{i}^{−1} mod q_{i} can be found with complexity GD(deg q_{i}) by the application to the pair Q_{i} mod q_{i}, q_{i} of the extended Euclidean algorithm which computes the linear representetion 1 = (Q_{i}, q_{i}) = UQ_{i}+V q_{i}, where U = Q_{i}^{−1} mod q_{i}, and p_{i} = (P mod q_{i})(Q_{i}^{−1} mod q_{i}) can be found with complexity 7 M(deg q_{i}) + 5 deg q_{i} by means of the multiplication and the subsequent division of the result by q_{i} (see [21, problems 17.32, 17.34]). With due account of superadditivity, we have the estimate GD(n) + 7 M(n) + 5n.

In order to calculate P mod q_{i}, we use the method of division in half. We construct a sequence of polynomials

Q_{α_{1},...,α_{k}}, α_{i} = 0, 1, k ≤ l(s) = ⌈log_{2} s⌉,

such that

Q = Q_{0} × Q_{1}, …, Q_{α_{1},...,α_{k}} = Q_{α_{1},...,α_{k},0} × Q_{α_{1},...,α_{k},1}, …,

and such that the number of factors in adjacent polynomials Q_{α̃,0}, Q_{α̃,1} would differ not more than by 1. Let us calculate the sequence of polynomials Pα̃ recursively:

P_{α̃,α_{k+1}} = Pα̃ mod Q_{α̃,α_{k+1}}.

If we set deg Qα̃ = nα̃, then the complexity of calculation of both polynomials P_{α̃,0}, P_{α̃,1} can be estimated as

3(M(n_{α̃,1}) + n_{α̃,1}) + 2 M(nα̃) + n_{α̃,0} + 3(M(n_{α̃,0}) + n_{α̃,0}) + 2 M(nα̃) + n_{α̃,1}
≤ 3 M(nα̃) + 4nα̃ + 4 M(nα̃) = 7 M(nα̃) + 4nα̃),

and the complete complexity can be estimated as

∑_{k=1}^{l(s)} ∑_{α_{1}=0,1;...;α_{k}=0,1} 7 M(nα̃) + 4nα̃ ≤ l(s)(7 M(  ∑_{α_{1}=0,1;...;α_{k}=0,1} nα̃) + 4 ∑_{α_{1}=0,1;...;α_{k}=0,1} nα̃)
≤ l(s)(7 M(n) + 4n).

With the use of Karatsuba’s method we evaluate the complexity of calculation of the sequence of polynomials Qα̃ as (l(s) − 1) M(n).

Since

Q_{i} mod q_{i} = \frac{Q_{i} mod q_{i}^{2}}{q_{i}} ,

this sequence can also be evaluated by means of dichotomy by the quantity of the same order of magnitude.

7. A more precise definition of multiplicative constants in the estimate of the complexity of the algorithm for factorization into partial fractions. Note that the sequence of remainders modulo Q2 ; α̃ that can arise in calculations is represented as Qα̃ ×Rα̃, where the sequence Rα̃ is computed recursively:

R_{α̃,α_{k+1}} = Q_{α̃,1−α_{k+1}} × Rα̃ mod Q_{α̃,α_{k+1}},

and the final results of Q_{i} mod q_{i} have the form R_{α_{1},...,α_{l(s)}}. Since deg Rα̃ < nα̃, the complexity of calculating both polynomials R_{α̃,0}, R_{α̃,1} is estimated (with the aid of problems 17.20, 17.25, 17.32 from [21]) as

M(nα̃) + 3(M(2n_{α̃,0}) + 2n_{α̃,0}) + M(nα̃ + n_{α̃,0}) + n_{α̃,1} + M(2n_{α̃,0}, n_{α̃,1}) + M(nα̃)
+ 3(M(2n_{α̃,1}) + 2n_{α̃,1}) + M(nα̃ + n_{α̃,1}) + n_{α̃,0} + M(2n_{α̃,1}, n_{α̃,0})
≤ 4 M(2nα̃) + M(3nα̃) + 3 M(nα̃) + 7nα̃,

and the complete complexity as l(s)(4 M(2n) + M(3n) + 3 M(n) + 7n). If we calculate the sequence Rα̃ before we calculate the sequence Pα̃, then, when calculating the latter sequence, we need not “invert” the polynomials Qα̃ since we have already performed this inversion when we calculated the preceding sequence (this is clear from the solution of problem 17.26 in [21]), and this makes it possible to lower the estimate of the complexity by 3 M(n) + 3n, and the total estimate will have the form l(s)(4 M(2n) + M(3n) + 8 M(n) + 8n) − M(n). We can lower the estimate still more if, before carrying out the multiplication

R_{α̃,α_{k+1}} = Q_{α̃,1−α_{k+1}} × Rα̃ mod Q_{α̃,α_{k+1}}

we reduce Rα̃ modulo Q_{α̃,α_{k+1}} and then “invert” the polynomial Qα̃ only once. Then we shall spend 2 M(n) operations for multiplication and 3 M(n) + 3n + 4 M(n) + 2n operations for modulo reduction. Then the total estimate will have the form l(s)(14 M(n) + 6n) − M(n). We can also note that, when the condition n_{α̃,0} = n_{α̃,1} of the equality of the degrees of the adjacent polynomials is fulfilled, this estimate decreases by 2l(s) M(n).

8. Estimate of the complexity of the second subalgorihm of the algorithm for decomposition of a rational fraction into partial fractions. Here we speak about the decomposition of a proper fraction of the form P/Q^{s} into the sum of still more simple fractions

∑_{i=1}^{s} \frac{P_{i}}{Q^{i}} .

Obviously, this problem is equivalent to the decomposition of the polynomial P in the powers of Q. It can be solved in the same way as the above-given computation of the remainders of a polynomial from the given moduli, only here all q_{i} are equal to Q and at every step the polynomial Pα̃ is used to construct the polynomials P_{α̃,0}, P_{α̃,1} such that

Pα̃ = P_{α̃,0}Q_{α̃,1} + P_{α̃,1}.

If, for every division, we choose the polynomial with the highest degree from two adjacent polynomials Qα̃, then the complexity of all divisions will be estimated (with the aid of problems 17.20, 17.25, 17.31, 17.32 from [21]) as

∑_{k=1}^{l(s)} ∑_{α_{1}=0,1;...;α_{k}=0,1}(3 M(n_{α̃,0}) + 3n_{α̃,0} + M(n_{α̃,1}) + n_{α̃,1} + 2 M(\frac{nα̃}{2} ) + \frac{nα̃}{2} )

≤ l(s)(4 M(\frac{n}{2}) + M(n) + \frac{5n}{2}),

and the complexity of calculation os all polynomials Qα̃, which result from a successive squaring or multiplication of adjacent polynomials, is evaluated as 2 M(n), where n = s deg Q, i.e., the complexity of the whole subalgorithm does not exceed l(s)(3 M(n) + 5n/2) + 2 M(n).

We see that Karatsuba’s method was actually used here as well.

The final estimate of the complexity of the algorithm for decomposition into partial fractions has the form

⌈log_{2} k⌉(17 M(n) + \frac{17n}{2}) + 8 M(n) + 5n + GD(n),

where k is the length of the decomposition of the denominator into squre-free factors, q_{1}q_{2}^{2} ⋯ q_{k}^{k}, and n = deg Q.

9. Estimate of the complexity of the algorithm proper for integration of rational fractions. We shall first estimate the complexity of integration of the fraction

\frac{P}{Q^{s}} = ∑_{i=1}^{s} \frac{P_{i}}{Q^{i}} ,

where deg P_{i} < deg Q, s deg Q = n, s > 1. For this purpose, we use the identity

∫ \frac{P_{i}}{Q^{i}} dx = \frac{−P_{i}B/(i − 1)}{Q^{i−1}} + ∫ \frac{P_{i}A + d(P_{i}B/(i − 1))/dx}{Q^{i−1}} dx,

where AQ+B ×(dQ/dx) = 1, deg A, deg B < deg Q, which is obtained by integration by parts. We evaluate the complexity of computing the polynomials A and B as GD(deg Q) and then find the numerator of the fraction under the sign of the second integral with complexity 2 M(deg Q)+6 deg Q, and divide it by Q with a remainder, representing it as CQ + D, where deg C, deg D < deg Q, with complexity 6 M(deg Q) + 5 deg Q, and then sum up C with P_{i−2} and D with P_{i−1} with complexity 2 deg Q. As a result, the problem reduces to a similar problem but i decreases by 1. This means that the integration of the fraction

\frac{P}{Q^{s}} = ∑_{i=1}^{s} \frac{P_{i}}{Q^{i}}

is carried out with complexity

(s − 1)(5 M(deg Q) + 10 deg Q) + 3 M(deg Q) + deg Q + GD(deg Q)

since we must compute A, B and “invert” the polynomial Q only once and need not divide by i − 1 when i = 2.

We must also have in mind the complexity of reduction of the rational part of the integral to the form

S_{0} + ∑_{i=1}^{s−1} \frac{S_{i}}{Q^{i}} ,

where deg S_{i} < deg Q. It should be noted that −P_{i}B/(i − 1) is computed with complexity M(deg Q)+2 deg Q (for i = 2 the last term vanishes), the division of the result by Q and, hence, its representation in the form C_{i}Q + D_{i}, where deg C_{i}, deg D_{i} < deg Q, is carried out with complexity 3 M(deg Q) + 2 deg Q (since “inversion” of the polynomial Q was performed earlier), and the summation (for i > 2) of the polynomials C_{i} and D_{i−1} is performed with complexity deg Q. The final result is the estimate

(s − 1)(4 M(deg Q) + 5 deg Q) − 3 deg Q.

Adding this estimate to that obtained earlier, we finally have the estimate

(s − 1)(9 M(deg Q) + 15 deg Q) + 3 M(deg Q) − 2 deg Q + GD(deg Q).

Using the superadditivity of the function M(n), we reduce this estimate to the form

6 M(n − deg Q) + 15n + 3 M(n) − 17 deg Q + GD(deg Q).

Applying this estiamte to each of the fractions resulting from the decomposition of the initial fraction N/P, where P = q_{1}q_{2}^{2} ⋯ q_{k}^{k}, and n = deg P, and using the superadditivity of the functions M(n) and GD(n), we get the estimate of the complexity of the integration proper in the form

3 M(n−deg q_{1})+6 M(p)+15p−2 deg Q+2 deg q_{1}+GD(deg Q) < 3 M(n)+6 M(p)+GD(n−p)+15p,

but this time we understand Q as the polynomial q_{1}q_{2} ⋯ q_{k}, and p = deg P_{1}, where P_{1} = P/Q.

10. Estimate of the complexity of addition of the obtained rational fractions. Ostrogradskii himself represented the result of integration as

∫ \frac{N}{P} dx = \frac{X}{P} + ∫ \frac{Y}{Q} dx.

In order to find the polynomials X and Y, we must sum up the obtained fractions (and the polynomials in the rational parts of the integrals). Applying the method of division in half and noting that the addition of two fractions of degree d requires not more than 3 M(d)+2d operations, we can compute the sum of proper fractions

∑_{i=1}^{s} \frac{P_{i}}{Q^{i}}

(solving the inverse of the problem considered at the end of Section 8) in

l(s)(M(s deg Q) + \frac{1}{2}s deg Q)

operations. Polynomials of a degree lower than deg q_{i} are added with complexity deg Q and then their sum is added to the proper fraction obtained earlier with complexity 2 deg Q + M(deg Q). Thus, for the complexity of addition of fractions we get the estimate

3 deg Q + M(deg Q) + ⌈log_{2} k⌉(3 M(p) + 2p + M(p) + \frac{p}{2} + 3 M(deg Q) + 2 deg Q)

≤ ⌈log_{k} t⌉(3 M(n) + M(p) + \frac{5p}{2}) + M(n − p) + 3(n − p).

The resulting estimate of the complexity of the algorithm for integration of a rational fraction has the form

⌈log_{2} k⌉(3 M(n) + M(p) + \frac{5p}{2}) + M(n − p) + 3(n − p) + 3 M(n) + 6 M(p) + GD(n − p) + 15p

+ ⌈log_{2} k⌉(17 M(n) + \frac{17n}{2}) + 8 M(n) + 5n + GD(n) + 2 GD(n) + 7 M(n) + 6n

≤ ⌈log_{k}(20 M(n) + M(p) + \frac{17n}{2} + \frac{5p}{2})

+ 3 GD(n) + GD(n − p) + 19 M(n) + 5 M(p) + 14n + 12p

≤ ⌈log_{2} n⌉(21 M(n) + 11n) + 24 M(n) + 26n + 4 GD(n).

11. Estimate of the complexity of finding the logarithmic part of an integral. The logarithmic part of the integral

∫ \frac{N}{P} dx,

namely,

∫ \frac{Y}{Q} dx,

cannot be found by purely arithmetic means (this problem is discussed, e.g., in [23]). However, if we know the roots of the polynomial Q, then it is, obviously, sufficient to decompose the fraction Y/Q into partial fractions in order to find the algorithmic part. Indeed, if Q = q_{1} ⋯ q_{m}, where q_{i}(x) = x − c_{i}, c_{i} ∈ ℂ, (q_{i}, q_{j}) = 1, i ≠ j, then we set Q_{i} = Q/q_{i} and, as is well known, obtain the decomposition

\frac{Y}{Q} = ∑_{i=1}^{s} \frac{p_{i}}{q_{i}} ,

where

p_{i} = Y × Q_{i}^{−1} mod q_{i} = Y mod q_{i}Q_{i}^{−1} mod q_{i} = \frac{Y (c_{i})}{Q′(c_{i})} = b_{i} ∈ ℂ.

Since the numbers Y (c_{i}) and Q′(c_{i}) are the remainders of the division of the plynomials Y (x) and Q′(x) = dQ/dx by the binomials q_{i}(x) = x − c_{i}, it follows, as was already noted in a more general situation, that they can be computed with complexity l(m)(15 M(m) + 8m) + 2m − M(m) (we have taken into account here that it suffices to calculate only once the products of the binomials q_{i} which are necessary for the computation of both numerators and denominators of the fractions b_{i}). Then we carry out the pairwise addition of the complex conjugate fractions with complexity not higher than 5m/2 and integrate all fractions in a simple way with complexity not exceeding 4m.

Note that if we are given the decomposition of the denominator P of the fraction being integrated into factors over the field ℝ, then there is no need to begin with calculating the polynomial Q by Ostsgradskii’s method since the problem reduces at once to the decomposition of the fraction N/P into the sum of fractions of the form Y/Q^{s}, where Q is either a linear binomial or a quadratic trinomial. In the first case, the decomposition of such a fraction into the sum of partial fractions reduces to finding the Taylor expansion of the polynomial P relative to the point a, where Q(x) = x − a, and, as is shown in [22], this can be done quicker than in the process described above, namely, with complexity O(s log s). The second case, as was noted in [17], reduces to the first one if we isolate a perfect square in Q, make the requisite linear change of variables in the polynomial Y (and this is the same as the Taylor expansion just mentioned), and then present it in the form Y_{0}(x^{2})+xY_{1}(x^{2}). After the substitution z = x^{2} the problem again reduces to the first case, and therefore, in the second case, the complexity of finding the decomposition of the fraction under consideration into the sum of partial fractions remains the same in the order of magnitude.

12. Estimate of the complexity of the original Ostrogradskii algorithm for integration of rational fractions. The algorithm of square-free factorization of polynomials presented above is only the beginning of Ostrogradskii’s algorithm for integration of rational fractions [1]. This algorithm is too large to be given here (the work [1] is 50 pages long). Using the techniques described above, we can evaluate its complexity as

2 GD(n) + GD(p) + O(\frac{nk M(n_{2})}{n_{2}}) + O(log k × M(n)),

where p = n − n_{1} and n_{i} = deg Q_{i}, and the sequence of polynomials Q_{i} is the same as that given in Section 3. Generally speaking, this estimate is worse than the estimate of the complexity of Hermite’s algorithm, but for p/n → 0 and (k log n_{2}) = o(log^{2} n) it is asymptotically twice as good.

13. Estimate of the complexity of the algorithm of the Chinese remainder theorem and of the problem of interpolation of polynomials. The Chinese remainder theorem and the related algorithms can boast of numerous applications in different divisions of mathematics and in many specific works (say, the algorithm for fast multiplicaton [8] uses modular arithmetic which is really based on the Chinese remainder theorem). It is interesting to note that the fast version of the Chinese algorithm is closely connected with the fast algorithm of decomposition into partial fractions, and both of them are actually based on the idea which was used for the first time by Karatsuba in [4].

Suppose that we are given polynomials q_{1}, …, q_{k}, (q_{i}, q_{j}) = 1, i ≠ j, and polynomials p_{1}, …, p_{k}, deg p_{i} < n_{i} = deg q_{i}, n = n_{1} + ⋯ + n_{k}. The Chinese remainder theorem asserts the existence and uniqueness of a polynomial P modulo Q such that P mod q_{i} = p_{i}, 1 ≤ i ≤ k. We set

Q = q_{1} ⋯ q_{k},  Q_{i} = \frac{Q}{q_{i}},  R_{i} = Q_{i}^{−1} mod q_{i},

and then, as is known, the relations

P = ∑_{i=1}^{s}(p_{i}Q_{i}R_{i}) mod Q = Q ∑_{i=1}^{s} \frac{(p_{i}R_{i}) mod q_{i}}{q_{i}}

holds true. Therefore, in order to compute the polynomial P, it suffices to find the products r_{i} = p_{i}R_{i} mod q_{i}, and this can be done in the same way as in Sections 6, 7 with complexity GD(n) + 7 M(n) + 5n + l(k)(10 M(n) + 5n) − M(n), and then compute the sum of proper fractions of degree n_{i}

\frac{P}{Q} = ∑_{i=1}^{k} \frac{r_{i}}{q_{i}}

just as in Section 10 with complexity l(k)(3 M(n) + n). The final estimate has the form

GD(n) + l(k)(13 M(n) + 6n) + 6 M(n) + 5n.

When q_{i}(x) = x − c_{i} and c_{i} and p_{i} belong to the field of coefficients, the polynomial P, which is constructed with the aid of the inverse Chinese algorithm, is known to be the solution of the interpolation problem P (c_{i}) = p_{i}. In this case, the computation of R_{i} = Q_{i}^{−1} mod q_{i} from the given Q_{i} mod q_{i} instead of the application of the extended algorithm for finding the GCD of polynomials requires only n division operations. It follows from what was said above that the complexity of its solution can be estimated as

l(n)(13 M(n) + 6n) + 6 M(n) + 6n.

If n is the power of 2, then, with due account of the remark made in Section 6, we can lower slightly this estimate.

The direct Chinese algorithm uses the polynomial P to compute its remainders from the given moduli q_{i} and, in particular, its values at the given points c_{i}. The estimate of its complexity was actually obtained in Section 6.

Note that all the results given above can be generalized to the numerical version of the Chinese remainder theorem and to the numerical version of the problem of the decomposition of a fraction into the sum of partial fractions. It should also be pointed out that in [13] the question concerning the fast Chinese algorithm (in the version without preliminary information processing) is presented not very clearly and the problem of the decomposition of a fraction into partial fractions is not considered at all (articles [17, 18] appeared later).

14. Estimate of the complexity of multiplication of polynomials and the algorithm of the Chinese remainder theorem. In order to compute the polynomial r, which is the product of the polynomials p and q of a degree smaller than n/2, it is natural to perform their multiplication modulo x^{n} − 1, for which purpose we must first calculate the remainders of the division of these polynomials by linear pairwise coprime binomials x − ε^{k} resulting from the factorization of the polynomial x^{n} − 1 (say, over the field ℂ), i.e., simply calculate the values p(ε^{k}) and q(ε^{k}), k = 0, …, n − 1, and then, multiplying them pairwise, obtain the values r(ε^{k}) = p(ε^{k})q(ε^{k}) and, using the “inverse” Chinese algorithm, restore the polynomial r. It is convenient to assume n to be equal to the power of 2 (adding, when necessary, the zero lowest coefficients to the polynomials p and q).

Applying the technique presented in Section 6 to compute the values p(ε^{k}) and q(ε^{k}) for k = 0, …, n − 1, we first find

p mod x^{n/2} − 1, p mod x^{n/2} + 1, q mod x^{n/2} − 1, q mod x^{n/2} + 1,

then calculate the remainders moduli x^{n/4} + 1, x^{n/4} − 1, x^{n/4} + i, x^{n/4} − i, and so on, until we find remainders moduli x^{2} − ε^{2k}, k = 0, …, n/2, x − ε^{k}, k = 0, …, n. Since the division of a polynomial of a degree lower than m by a binomial of degree m/2 is carried out by a “school” algorithm with complexity m, the complexity of the whole algorithm for computing the values p(ε^{k}) and q(ε^{k}), k = 0, …, n − 1, is estimated as 2n log_{2} n (multiplications and divisions carried out in the field ℂ).

For the restoration (interpolation) of the polynomial r from its known values r(ε^{k}), k = 0, …, n − 1, it is natural to use the Lagrange formula

r = f(x) ∑_{k=0}^{n−1} \frac{r(ε^{k})/f′(ε^{k})}{x − ε^{k}} = f(x) ∑_{k=0}^{n−1} \frac{r(ε^{k})ε^{k}}{n(x − ε^{k})},

where f(x) = x^{n} − 1, and perform the summation of the fractions in the same way as we did in Section 10. Since, when the order of summation is properly chosen, all the resulting fractons have binomial denominators and the multiplication of a polynomial of a degree lower than m by the “school” method is of complexity m, the addition of two fractions of this kind has the complexity 4m, and the whole interpolation algorithm has the complexity 2n + 2n log_{2} n. As a result, we get the well-known relation M(n) = O(n log n).

We see that the fast multiplication of polynomials is actually carried out with the aid of the Chinese remainder theorem and the method of division in half which was first applied to this problem by Karatsuba [4]. As Karatsuba recounts in his articale [5], Kolmogorov spoke of the possibility of application of the Chinese remainder theorem to multiplication of numbers and initiated the research in this field.

The descrete Fourier transformation was not mentioned in the given proof, although, actually, the caculation of the values of a polynomial in nth degree roots of unity, which was carried out above, coincides with the discrete Fourier transformation, and the technique, which can be used for this purpose, coincides with one of the variety of algorithms of Fourier fast transformations known today (the monograph [24] is devoted to these algorithms). A similar algorithm for Fourier fast transformation is described in [13]. The interpolation of a polynomial from the values in nthdegree roots of unity coincides with Fourier inverse transformation, and therefore, the algorithm indicated above could not be used for its performance since the estimate of the complexity of the inverse Fourier transformation coincides, as is well known, with the estimate of the copmplexity of the direct Fourier transformation. The application of the Chinese theorem in this situation is equivalent to the application of the so-called theorem on a cyclic convolution.

ACKNOWLEDGMENTS

This work was supported by the Russian Foundation for Fundamental Research (project no. 97-01-00721).

REFERENCES

1. Ostrogradskii, M. V., On the Integration of Rational Fractions, in: Selected Works, Moscow: Akad. Nauk SSSR, 1958, pp. 56–105.

2. Ostrogradskii, M. V., Note on the Equal Factors of Integral Polynomials, in: Selected Works, Moscow: Akad. Nauk SSSR, 1958, pp. 106–107.

3. Hermite, C., Quevres de Charles Hermite, vol. 3 (Picard, E., Ed.), Paris: Gautnier-Villars, Imprimeur-Librae, 1912.

4. Karatsuba, A. A. and Ofman, Yu. P., Multiplication of Multiplace Numbers with the Use of Automata, Dokl. Akad. Nauk SSSR, 1962, vol. 145, no. 2, pp. 293–294.

5. Karatsuba, A. A., Complexity of Computations, Trudy Mat. Inst. Rus. Akad. Nauk, 1995, vol. 211, pp. 186–202.

6. Toom, A. L., On the Complexity of the Scheme of Functional Elements Which Realizes the Multiplication of Integers, Dokl. Akad. Nauk SSSR, 1963, vol. 150, no. 2, pp. 496–498.

7. Schönhage, A. and Strassen, V., Schnelle Multiplikation grossen Zahlen, Computing, 1971, vol. 7, no. 3–4, pp. 281– 292.

8. Schönhage, A., Schnelle Multiplikation grossen Zahlen, Computing, 1966, vol. 1, pp. 182–196.

9. Cook, S. A., On the Minimum Computation Time of Functions, Thesis Harvard Univ., 1966, pp. 51–77.

10. Benderskii, Yu. V., Fast Computations, Dokl. Akad. Nauk SSSR, 1975, vol. 223, no. 5, pp. 1041–1043.

11. Kung, H. T., Fast Evaluation and Interpolation, Pittsburgh: Dept. Comput. Sci., Carnegie-Mellon Univ., 1973.

12. Strassen, V., Gaussian Elimination is not Optimal, Numer. Math., 1969, vol. 13, no. 4, pp. 354–356.

13. Aho, A., Hopcroft, J., and Ullman, J., The Design and Analysis of Computer Algorithms, USAM: Addison-Wesley, 1975.

14. Knuth, D., The Art of Computer Programming, vol. 2, Reading: Addison-Wesley, 1969.

15. Schönhage, A., Grotefeld, A. F. W., and Vetter, E., Fast Algorithms, B. 1, Manheim, Leipzig, Wien, Zürich: Wissenschaftsverlag, 1994.

16. Yun, D. Y. Y., Fast Algorithm for Rational Function Integration, in: Proc. IFIP Congr., Amsterdam: North-Holland, 1977, vol. 7, pp. 493–498.

17. Kung, H. T. and Tong, D. M., Fast Algorithms for Fraction Decompositon, SIAM J. Comput., 1977, vol. 6, pp. 582–593.

18. Chin, P. I., Complexity of Partial Fraction Expansion, SIAM J. Comput., 1977, vol. 6, pp. 554–562.

19. Moenk, R., Fast Algorithm of GCD’s, in: Proc. 5th Ann. ACM Sympos. on Theory of Computing, 1973, pp. 142– 151.

20. Schönhage, A., Schnelle berechnung von kettenbruchentwicklungen, Acta Inform., 1971, vol. 1, pp. 139–144.

21. Gashkov, C. B. and Chubarikov, V. N., Arifmetika, Algoritmy. Slozhnost Vychislenii (Arithmetic Algorithms. Complexity of Computations), Moscow: Nauka, 1996.

22. Aho, A. V., Steiglitz, K., and Ullman, J. D., Evaluation of Polynomials at Fixed Set of Points, SIAM J. Comput., 1975, vol. 6, pp. 533–539.

23. Davenport, J., On the Integration of Algebraic Functions, Heidelberg: Springer, 1981.

24. Blahut, R., Bystrye algoritmy tsifrovoi obrabotki signalov (Fast Algorithms of Digital Signal Proc.), Moscow: Mir, 1989.

Translated by I. Aleksandrova