Power series inversion and division

A power series Q = Q_0 + Q_1 x + \cdots over \mathbb{Z}/n\mathbb{Z} is invertible if and only if Q_0 is a unit, and then its inverse modulo x^m is a polynomial of length m. These functions compute 1/Q and A/B modulo x^m. Inversion is the basic step of division with remainder (applied to the reversed divisor, see Computing modulo a fixed polynomial) and of most power series functions, and it costs only a small constant number of multiplications.

Algorithms#

The basecase algorithm solves Q \cdot \mathit{Qinv} = 1 coefficient by coefficient:

with each sum computed as one dot product (_nmod_vec_dot_rev, see Computing dot products) and a single reduction. It costs O(m \min(m, \mathit{Qlen})) operations. A binomial Q_0 + c x^k is inverted as a geometric series. Division uses the analogous recurrence q_i = B_0^{-1} (A_i - \sum_{j \ge 1} B_j q_{i-j}).

Newton iteration doubles the number of correct coefficients in each step: if g \equiv 1/Q \pmod{x^k}, then

The product Qg is known to be 1 in its low k coefficients, so only its middle coefficients are needed, which is a middle product, and the correction is a truncated product of length k. Summed over all steps the cost is a small multiple of one multiplication of length m. Division A/B is computed in a similar way, folding the multiplication by A into the last Newton step instead of computing 1/B first, as proposed by Karp and Markstein [KM1997].

nmod_poly_inv_series uses the basecase algorithm for \mathit{Qlen} \le 10, and otherwise for m below a cutoff tabulated by the bit size of n (with fft_small: 36 to 50 for moduli up to 7 bits, 81 to 292 for 8 to 21 bits, 388 to 407 for 22 to 28 bits and 231 to 448 above) or for short Q, and Newton iteration otherwise. nmod_poly_div_series uses the basecase for \mathit{Blen} \le 5 and analogous tables (about 80 to 600). nmod_poly_inv_series_newton is the same function as nmod_poly_inv_series.

Examples#

The generating function of the Fibonacci numbers is x / (1 - x - x^2). One series division gives F_0, \ldots, F_{m-1} modulo p, and here F_{100000} \bmod p for p = 10^9 + 7, checked against the matrix-free identity F_{2k} = F_k (2 F_{k+1} - F_k):

#include <flint/flint.h>
#include <flint/nmod.h>
#include <flint/nmod_poly.h>

int main(void)
{
    nmod_poly_t num, den, f;
    slong m = 100001;
    ulong a, b;
    nmod_t mod;

    nmod_init(&mod, 1000000007);
    nmod_poly_init_mod(num, mod);
    nmod_poly_init_mod(den, mod);
    nmod_poly_init_mod(f, mod);

    nmod_poly_set_coeff_ui(num, 1, 1);                  /* x */
    nmod_poly_set_coeff_ui(den, 0, 1);                  /* 1 - x - x^2 */
    nmod_poly_set_coeff_ui(den, 1, mod.n - 1);
    nmod_poly_set_coeff_ui(den, 2, mod.n - 1);

    nmod_poly_div_series(f, num, den, m);

    flint_printf("F_0..F_10:");
    for (slong i = 0; i <= 10; i++)
        flint_printf(" %wu", nmod_poly_get_coeff_ui(f, i));
    flint_printf("\n");

    a = nmod_poly_get_coeff_ui(f, 50000);
    b = nmod_poly_get_coeff_ui(f, 50001);
    flint_printf("F_100000  = %wu\n", nmod_poly_get_coeff_ui(f, 100000));
    flint_printf("identity  = %wu\n",
        nmod_mul(a, nmod_sub(nmod_add(b, b, mod), a, mod), mod));

    nmod_poly_clear(num);
    nmod_poly_clear(den);
    nmod_poly_clear(f);
    return 0;
}

F_0..F_10: 0 1 1 2 3 5 8 13 21 34 55
F_100000  = 911435502
identity  = 911435502

The Bernoulli numbers are defined by x / (e^x - 1) = \sum_k B_k x^k / k!. Inverting the series (e^x - 1)/x = \sum_k x^k / (k+1)! modulo a prime p larger than the length gives B_k \bmod p; the program prints 6 B_2, 30 B_4, 42 B_6, 30 B_8 and 66 B_{10}, which are 1, -1, 1, -1, 5:

#include <flint/flint.h>
#include <flint/nmod.h>
#include <flint/nmod_poly.h>

int main(void)
{
    nmod_t mod;
    nmod_poly_t h, b;
    slong k, m = 12;
    ulong fac = 1, den[6] = { 0, 6, 30, 42, 30, 66 };

    nmod_init(&mod, 1000003);
    nmod_poly_init_mod(h, mod);
    nmod_poly_init_mod(b, mod);

    /* h = sum x^k / (k+1)! */
    for (k = 0; k < m; k++)
    {
        fac = nmod_mul(fac, k + 1, mod);
        nmod_poly_set_coeff_ui(h, k, nmod_inv(fac, mod));
    }

    nmod_poly_inv_series(b, h, m);

    /* B_k = k! * b_k */
    fac = 1;
    for (k = 1; k < m; k++)
    {
        fac = nmod_mul(fac, k, mod);
        if (k % 2 == 0)
        {
            ulong B = nmod_mul(fac, nmod_poly_get_coeff_ui(b, k), mod);
            ulong t = nmod_mul(B, den[k / 2], mod);
            flint_printf("%wu * B_%wd = %wd\n", den[k / 2], k,
                         (t > mod.n / 2) ? (slong) t - (slong) mod.n : (slong) t);
        }
    }

    nmod_poly_clear(h);
    nmod_poly_clear(b);
    return 0;
}

6 * B_2 = 1
30 * B_4 = -1
42 * B_6 = 1
30 * B_8 = -1
66 * B_10 = 5