Newton iteration for power series and division

Newton's method doubles the number of correct terms of a power series approximation with each step. Since the last step dominates, computing an inverse, a quotient, a square root or an exponential to n terms costs a small constant times one multiplication of length n. Power series inversion by Newton iteration is due to Sieveking [Sie1972] and Kung [Kun1974]; Brent and Kung [BrentKung1978] treated general equations and elementary functions. Textbook treatments are [GG2003, Ch. 9] and [BZ2011, §4.2]; Hanrot and Zimmermann [HZ2004] and Bernstein [Ber2008] discuss the constants.

Throughout, the coefficients lie in R = \mathbb{Z}/p\mathbb{Z} for a word-size modulus p (we write p rather than n to keep n for the precision), \beta is the bit length of p, "to precision n" means modulo x^n, and M(n) is the cost of multiplying two polynomials of length n. We assume, as is true for all methods in Multiplying polynomials over Z/nZ, that M(n)/n is nondecreasing, so \sum_{i \ge 1} M(\lceil n/2^i \rceil) \le M(n) + O(\log n). With the FFT, a low product of length n and a middle product with output length n cost about M(n) each (Middle products and the transposition principle).

Generic version#

The iterations described here are implemented once in gr_poly and are valid over any commutative ring, not only \mathbb{Z}/p\mathbb{Z}: each step uses ring operations and one inversion of the constant coefficient (the leading coefficient of the divisor for division with remainder), and the algorithms return GR_DOMAIN if that coefficient is not invertible. Over real and complex balls they give enclosures of the exact result, since every operation does. The generic dispatchers use the basecase for short inputs or when the ring has only classical multiplication (inverse and quotient of power series below length 9, division when the quotient or the divisor has length at most 10, exponential below length 10), Newton iteration otherwise, and middle products when the ring provides them (Division and Newton iteration). When a leading coefficient is not invertible, as for division by 2x + 1 over \mathbb{Z}, division with remainder and power series division fall back to algorithms that divide exactly at each step and succeed whenever the quotient lies in the ring. The rings nmod, fmpz_mod, mpn_mod and the fq family override the dispatchers to use their own basecase kernels and cutoffs, as nmod does with the tables quoted below; the ball rings divide through power series with error bounds.

The general step#

FLINT runs Newton iterations from the target precision downwards: it records a_0 = n, a_{i+1} = \lceil a_i / 2 \rceil until a_k is below a cutoff, computes the solution to precision a_k by a basecase (quadratic) algorithm, and then lifts from a_{i+1} to a_i. This avoids overshooting: no step computes more terms than are needed. The same scheme lifts p-adic solutions, with p in place of x (Hensel lifting of inverses, square roots and Teichmüller lifts).

Inverse#

Let Q(0) be a unit and g \equiv 1/Q \pmod{x^m}. Then e = 1 - Qg \equiv 0 \pmod{x^m}, and since Qg = 1 - e,

so g(1 + e) is correct to 2m terms and g(1 + e + e^2) to 3m. The first is the Newton step of Lemma 1 for \Phi = 1/y - Q. The low m coefficients of Qg are 1, 0, \ldots, 0, so only the window w = (Qg)_{[m, n)} has to be computed, which is a middle product, and e = -x^m w.

Exploiting the structure of the FFT further (reusing transforms between steps and between the two products) lowers the constant; the best known value is \tfrac{13}{9}M(n), due to Harvey [Har2011], who also gives \tfrac{4}{3}M(n) for the square root. FLINT does not use these FFT-specific methods.

The third-order step costs one more short squaring but triples the precision, so there are fewer steps. FLINT uses it for n \ge 8192 in rings of finite characteristic, where it is 5–10% faster with the FFT. (Over \mathbb{Z} and \mathbb{Q} the coefficients of e^2 can be much larger than those of the result.)

In FLINT. nmod_poly_inv_series uses the recurrence (with _nmod_vec_dot_rev and delayed reduction) for \operatorname{len}(Q) \le 10, and otherwise calls the gr_poly code with the nmod context, which uses the recurrence if n is below a cutoff depending on \beta (NMOD_BITS(mod)), or if Q is short, and _gr_poly_inv_series_newton otherwise. With fft_small the cutoff is about 40 for \beta \le 7, 80–300 for 8 \le \beta \le 21, about 400 for 22 \le \beta \le 28, and 230–450 above. The recurrence is used for any n when \operatorname{len}(Q) is below a second table (between 39 and 141 for \beta \ge 23; smaller moduli do not use this rule), since it costs only O(n \operatorname{len}(Q)). If Q is a binomial Q_0 + c x^k, the inverse is a geometric series and is written down directly. nmod_poly_inv_series_newton is the same function.

Quotient of power series#

To compute A/B \bmod x^n, one could invert B to n terms and multiply. The Karp–Markstein trick [KM1997] instead merges the multiplication by A into the last Newton step.

This costs an inverse to n/2 terms plus three products of length about n/2, about M(n) + \tfrac{3}{2}M(n) = \tfrac{5}{2}M(n), instead of 2M(n) + M(n) = 3M(n).

In FLINT. nmod_poly_div_series uses the recurrence q_i = B_0^{-1}(A_i - \sum_{j \ge 1} B_j q_{i-j}) for \operatorname{len}(B) \le 5, and otherwise _gr_poly_div_series_newton (Algorithm 2) above a cutoff that, with fft_small, is about 80–340 for \beta \le 20, about 600 for 21 \le \beta \le 28 and 320–490 above, with a separate table for short B.

Division with remainder#

Let A = QB + R with \deg A = a, \deg B = b, \deg R < b and \deg Q = a - b. Write \operatorname{rev}_k(f) = x^k f(1/x).

So the quotient is a power series quotient of length \ell_Q = a - b + 1, which needs only the top \ell_Q coefficients of A and of B, and the remainder is R = (A - QB) \bmod x^b, a low product.

By Lemma 4 the algorithm is correct. For \ell_A = 2n and \ell_B = n it costs about \tfrac{5}{2}M(n) + M(n) = \tfrac{7}{2}M(n), compared with n^2 operations for schoolbook division.

If many dividends are reduced modulo the same B, the power series inverse of \operatorname{rev}(B) can be computed once (it is usually computed to \ell_B terms; \ell_Q terms suffice). Then each quotient with \ell_Q \le \ell_B - 1 is one low product, and the division costs two products: this is how _nmod_poly_div_newton_n_preinv, _nmod_poly_divrem_newton_n_preinv and everything built on them (nmod_poly_mulmod_preinv, nmod_poly_powmod_ui_binexp_preinv, nmod_poly_compose_mod_brent_kung_preinv) work. For a sparse divisor (at most five nonzero coefficients below the leading one, and \ell_B > 20), _nmod_poly_divrem_newton_n_preinv instead eliminates directly, at a cost of a few operations per coefficient of the quotient.

In FLINT. nmod_poly_divrem and nmod_poly_rem use schoolbook division (_nmod_poly_divrem_basecase_preinv1, with delayed reduction) if \ell_A - \ell_B < 80 or \ell_B is less than 400 (for \beta \le 27 on x86-64 with AVX2; 200 without AVX2) or 80 (for \beta > 27), and Algorithm 3 otherwise. Without fft_small, divide-and-conquer division is used instead for \beta \ge 16, \ell_B \ge 1024 and \ell_A \le 16384. nmod_poly_div uses schoolbook division if \ell_B \le 15 or \ell_A - \ell_B \le 15, and steps 1 and 2 of Algorithm 3 otherwise.

Square roots#

Let 2 be invertible and let h(0) have a square root in R. For the reciprocal square root, let g \equiv h^{-1/2} \pmod{x^m} and e = 1 - hg^2 = O(x^m). Since h = g^{-2}(1 - e),

FLINT uses this third-order iteration throughout (_gr_poly_rsqrt_series_newton): each step computes the square g^2 (low product), the window of h g^2 (middle product), the short square e^2 and one low product with g, and triples the precision. When n \le 2m the term in e^2 is dropped.

The square root itself is obtained from the reciprocal square root to half precision by a Karp–Markstein step:

For short h, FLINT uses instead J. C. P. Miller's recurrence for powers [Knu1997, §4.7]: y = h^\alpha satisfies h y' = \alpha h' y, which gives each coefficient of y from the previous ones in O(\operatorname{len}(h)) operations. It divides by 1, \ldots, n - 1, so it is used only when n is less than the modulus, and only for \operatorname{len}(h) below a tuned bound (with fft_small, about 40 to 190 depending on \beta, and only for \beta \ge 19) and 8 \operatorname{len}(h) < n (square root) or 2 \operatorname{len}(h) < n (reciprocal square root).

nmod_poly_sqrt_series and nmod_poly_invsqrt_series call these algorithms. nmod_poly_sqrt, the square root of a polynomial that is a perfect square, removes an even power of x, computes the series square root to half the length plus one, and checks the result by comparing the high half of its square (_nmod_poly_mulhigh) with the input.

Exponential and logarithm#

For h(0) = 0, \log and \exp of power series are defined if 1, 2, \ldots, n - 1 are invertible, which for prime modulus p means n \le p. The logarithm is

one power series division and an integration (_gr_poly_log_series, called by nmod_poly_log_series; if f is a binomial f_0 + c x^k the coefficients are written down directly).

For the exponential, FLINT uses a Newton iteration that computes f = \exp(h) and g = \exp(-h) = 1/f together [HZ2004; BZ2011, §4.2.5].

Each step costs two middle products and three low products of length about n/2; the last step skips the update of g unless \exp(-h) is requested. Faster but more complicated methods are known: \tfrac{7}{3}M(n) by van der Hoeven [vdH2010] and \tfrac{13}{6}M(n) by Harvey [Har2009b], both relying on FFT-specific techniques.

The basecase uses the recurrence k f_k = \sum_{j=1}^{k} j h_j f_{k-j} from f' = h' f, costing O(n \operatorname{len}(h)). An intermediate method (_gr_poly_exp_series_basecase_mul) splits h = a + x^m b with m = \lceil n/3 \rceil and uses \exp(h) \equiv \exp(a)(1 + x^m b + x^{2m} b^2/2) \pmod{x^n}. nmod_poly_exp_series uses the basecase if \operatorname{len}(h) is below a tuned bound (230–490 with fft_small), the split method below a second bound (470–1491), and Algorithm 4 otherwise, with the second bound as its basecase cutoff. nmod_poly_exp_expinv_series computes the exponential and then inverts it; nmod_poly_sinh_series and nmod_poly_cosh_series use it. The other trigonometric and inverse trigonometric series functions are computed by their gr_poly counterparts.

Use in FLINT#

All the Newton iterations live in gr_poly and are called by nmod_poly through the gr context for \mathbb{Z}/n\mathbb{Z}, which supplies the nmod_poly multiplication and middle product and the tuned cutoffs quoted above (tables indexed by \beta).

The nmod_poly cutoffs between the basecase algorithms and Newton iteration are listed in Algorithms for series inversion and division and in Algorithms for division with remainder.

Power series composition and reversion, which also use Newton iteration for long series, are treated in Composition and modular composition.