Hermite reduction

Hermite reduction computes the rational part of the integral of a rational function^[1]^[4]. For $A, D \in K [x]$ over a field $K$ of characteristic zero, the antiderivative of $A / D$ is a rational function plus a sum of constant multiples of logarithms, and the rational function comes entirely from the repeated factors of $D$ . Hermite's method extracts it using the square-free decomposition of $D$ and the extended Euclidean algorithm, never factoring $D$ into irreducibles, and leaves an integral with a square-free denominator. The logarithmic part of that remainder is then given by the Rothstein–Trager resultant^[5]^[6]. The same reduction, applied in the top variable of a tower of extensions, is the first step of the Risch algorithm. Ostrogradsky found an equivalent method in 1845^[2] and Hermite published his in 1872^[1].

§ 01

Example

\int \frac{d x}{( x ^{2} + 1 ) ^{2}} = \frac{x}{2 ( x ^{2} + 1 )} + \int \frac{1/2}{x ^{2} + 1} d x = \frac{x}{2 ( x ^{2} + 1 )} + \frac{1}{2} arctan x

The denominator has a repeated factor $x^{2} + 1$ . Hermite reduction produces the rational term $x / (2 (x^{2} + 1))$ and leaves an integrand whose denominator $x^{2} + 1$ is square-free. The remainder integrates to logarithms over $Q (i)$ : $arctan x = \frac{i}{2} (lo g (x + i) - lo g (x - i))$ .

§ 02

The reduction step

Divide first, so that $de g A < de g D$ ; the polynomial quotient integrates term by term. Write $D = U V^{m}$ with $V$ square-free, $g cd (U, V) = 1$ and $m \geq 2$ . Since $V$ is square-free, $g cd (V, V^{'}) = 1$ , so $g cd (U V^{'}, V) = 1$ and the extended Euclidean algorithm gives $B$ and $C$ with^[1]^[4]

U V^{'} B + V C = A

and then

\int \frac{A}{U V ^{m}} = \frac{B}{( 1 - m ) V ^{m - 1}} + \int \frac{C - U B ^{'} / ( 1 - m )}{U V ^{m - 1}}

To check it, differentiate the first term:

(\frac{B}{( 1 - m ) V ^{m - 1}})^{'} = \frac{B ^{'}}{( 1 - m ) V ^{m - 1}} + \frac{B V ^{'}}{V ^{m}}

Adding the integrand on the right, the $B^{'}$ terms cancel and what is left is $(U B V^{'} + C V) / (U V^{m}) = A / (U V^{m})$ . Each step lowers the power of $V$ by one. Taking $V = D_{k}$ from the square-free decomposition $D = D_{1} D_{2}^{2} \dots D_{k}^{k}$ and repeating ends with a square-free denominator.

For $\int d x / (x^{2} + 1)^{2}$ : $U = 1$ , $V = x^{2} + 1$ , $m = 2$ , $V^{'} = 2 x$ . Solving $2 x B + (x^{2} + 1) C = 1$ gives $B = - x /2$ and $C = 1$ , and the step produces $x / (2 (x^{2} + 1))$ plus the integral of $(1 - \frac{1}{2}) / (x^{2} + 1)$ , as above.

§ 03

The split it computes

Hermite reduction splits the square-free decomposition of D. One copy of each factor, the square-free part D*, stays under the integral; the other copies form D⁻ = gcd(D, D′), the denominator of the rational part.

The result of the reduction is determined in advance. With $D^{-} = g cd (D, D^{'})$ and $D^{*} = D / D^{-}$ , there are unique polynomials $P$ and $Q$ with $de g P < de g D^{-}$ and $de g Q < de g D^{*}$ such that^[2]^[3]

\int \frac{A}{D} = \frac{P}{D ^{-}} + \int \frac{Q}{D ^{*}}

In terms of the square-free decomposition, $D^{-} = D_{2} D_{3}^{2} \dots D_{k}^{k - 1}$ and $D^{*} = D_{1} D_{2} \dots D_{k}$ . Uniqueness holds because the integral of a fraction with square-free denominator has no rational part: a rational function whose derivative had only simple poles would have to be a polynomial.

Ostrogradsky and Horowitz compute $P$ and $Q$ directly^[2]^[3]. Differentiating the identity and clearing denominators gives

A = P^{'} D^{*} - P \frac{D ^{*} D ^{-'}}{D ^{-}} + Q D^{-}

which is linear in the unknown coefficients of $P$ and $Q$ , a system of $de g D$ equations in $de g D$ unknowns. Solving it by Gaussian elimination costs $O (n^{3})$ for $n = de g D$ , which is more than Hermite's sequence of Euclidean steps^[4]^[8].

§ 04

Algorithm

Hermite's original method needs the square-free decomposition and a partial fraction split across its factors. A variant attributed by Bronstein to Mack reduces all the factors of the highest remaining multiplicity at once, using only $D^{-}$ and $D^{*}$ and a gcd per step, and each Diophantine equation it solves has size proportional to the part being reduced^[4]. Writing $D_{m} = g cd (D, D^{'})$ and $D_{s} = D / D_{m}$ , each round takes $D_{m 2} = g cd (D_{m}, D_{m}^{'})$ and $D_{m s} = D_{m} / D_{m 2}$ , the product of the factors of highest remaining multiplicity, solves

B \cdot (- \frac{D _{s} D _{m}^{'}}{D _{m}}) + C \cdot D_{m s} = A

adds $B / D_{m}$ to the rational part, replaces $A$ by $C - B^{'} D_{s} / D_{m s}$ and $D_{m}$ by $D_{m 2}$ , and stops when $D_{m}$ is constant. The remaining integrand is $A / D_{s}$ with $D_{s} = D^{*}$ square-free.

Hermite reduction in this form, over the rationals.

(* Hermite reduction, in Mack's linear form (Bronstein, Algorithm 2.2.2):
   returns the polynomial part q, the rational part g and a fraction a/ds
   with integral(A/D) = integral(q) + g + integral(a/ds), where ds is
   square-free and deg a < deg ds. g is returned as numerator and
   denominator. *)
let hermite (a : t) (d : t) =
  let q, a = divmod a d in
  let dm = gcd d (deriv d) in
  let ds = quo d dm in
  let rec go a dm (gn, gd) =
    if deg dm <= 0 then ((gn, gd), a)
    else begin
      let dm2 = gcd dm (deriv dm) in
      let dms = quo dm dm2 in
      (* solve b (-ds dm' / dm) + c dms = a *)
      let b, c = diophantine (neg (quo (mul ds (deriv dm)) dm)) dms a in
      let a' = sub c (quo (mul (deriv b) ds) dms) in
      (* g <- g + b / dm *)
      let gn' = add (mul gn dm) (mul b gd) and gd' = mul gd dm in
      let h = gcd gn' gd' in
      go a' dm2 (quo gn' h, quo gd' h)
    end in
  let g, a = go a dm ([||], one) in
  (q, g, (a, ds))
 
let rec qpow c n = if n = 0 then Q.one else Q.mul c (qpow c (n - 1))

Two reductions, each checked by differentiating the rational part.

open Qpoly

Running it.

1 / (x^2 + 1)^2
  rational part  ((1/2)x)/(x^2 + 1)
  remaining      integral of (1/2)/(x^2 + 1)
  derivative checks: true
(x^4 + 1) / ((x + 1)^3 (x^2 + 1)^2)
  rational part  (-(1/4)x^3 - (1/2)x^2 + (1/4)x)/(x^4 + 2x^3 + 2x^2 + 2x + 1)
  remaining      integral of (-(1/4)x + 3/4)/(x^3 + x^2 + x + 1)
  derivative checks: true

Both agree with SymPy's ratint_ratpart. The second example has two repeated factors, of multiplicities 3 and 2; the remaining denominator $x^{3} + x^{2} + x + 1 = (x + 1) (x^{2} + 1)$ is their product.

§ 05

No factorization needed

The reduction uses only arithmetic, gcds and the extended Euclidean algorithm in $K [x]$ , so its output has coefficients in $K$ even when the factors of $D$ are irreducible over $K$ or have no roots in it^[4]. For a denominator like $(x^{3} + x + 1)^{3}$ , whose roots are the roots of an irreducible cubic, the rational part is computed without ever naming them.

The rational part of 1/(x³ + x + 1)³, then the Rothstein–Trager resultant of the remainder.

open Qpoly

Running it.

D = x^9 + 3x^7 + 3x^6 + 3x^5 + 6x^4 + 4x^3 + 3x^2 + 3x + 1
rational part  (-(243/961)x^5 + (357/1922)x^4 - (405/961)x^3 - (315/1922)x^2 + (156/961)x - 224/961)/(x^6 + 2x^4 + 2x^3 + x^2 + 2x + 1)
remaining      integral of (-(243/961)x + 357/961)/(x^3 + x + 1)
Rothstein-Trager R(z) = z^3 - (1104057/28629151)z - 84213/28629151
R has no rational roots: the logarithms need algebraic numbers

The rational part agrees with SymPy. The logarithmic part here needs the three roots of a cubic in $z$ , which are algebraic numbers of degree 3; the answer can be written as a sum over those roots without computing them.

§ 06

The logarithmic part

After reduction the integrand $A / D$ has $D$ square-free and $de g A < de g D$ . Over the algebraic closure it equals $\sum_{α} res_{α} / (x - α)$ with residues $A (α) / D^{'} (α)$ at the roots of $D$ , so its integral is $\sum_{α} res_{α} lo g (x - α)$ . Grouping the roots by residue gives the Rothstein–Trager theorem^[5]^[6]: with

R (z) = res_{x} (D, A - z D^{'})

the distinct roots $c_{i}$ of $R$ are exactly the distinct residues, and

\int \frac{A}{D} = i \sum c_{i} lo g g cd (D, A - c_{i} D^{'})

The logarithms are taken of the smallest possible arguments: the algebraic extension of $K$ needed is the one generated by the roots of $R$ , which is often much smaller than the splitting field of $D$ . Lazard and Rioboo, and independently Trager, showed that the gcds can be read off the subresultant sequence of $D$ and $A - z D^{'}$ , avoiding computation with algebraic numbers entirely^[7].

A full rational integrator: polynomial part, Hermite reduction and Rothstein–Trager, for cases where the residues are rational.

open Qpoly

Running it.

integral of (x^6 + 2x^4 - x^3 + 1)/(x^5 - 2x^4 + x^3)
  = (1/2)x^2 + 2x + (-5x^2 + (3/2)x + 1/2)/(x^3 - x^2) + 2 log(x - 1) + 3 log(x)
integral of (x^2 + 3)/(x^4 + 4x^3 + 3x^2 - 4x - 4)
  = (-7/3)/(x + 2) - 2 log(x + 1) + 2/9 log(x - 1) + 16/9 log(x + 2)

Both answers agree with SymPy's ratint.

§ 07

Complexity

For a denominator of degree $n$ , Hermite reduction in the form above uses $O (n^{2})$ operations in $K$ with classical arithmetic^[4]^[8]. Over $Q$ the coefficients of intermediate results grow, and Gerhard gave modular algorithms for square-free factorization and Hermite integration with nearly optimal running time in both the degree and the coefficient size^[9].

§ 08

History

Ostrogradsky described the split into $P / D^{-}$ and a square-free remainder in 1845, computing $P$ and $Q$ from a linear system^[2]. Hermite gave the reduction by Euclidean steps in 1872^[1]. Both methods returned in computer algebra: Horowitz analyzed the cost of both approaches for rational function integration in 1971^[3], and Rothstein and Trager independently found the resultant that gives the logarithmic part without factoring in 1976^[5]^[6]. Lazard and Rioboo gave the subresultant method for the logarithmic arguments in 1990^[7]. Bronstein's book presents the modern forms of all these algorithms and their generalization to transcendental extensions, where they form the rational-function steps of the Risch algorithm^[4].