En travaillant avec des polynômes, on constate rapidement que les complexités algorithmiques peuvent croître de manière exponentielle lors d'appels de fonctions imbriqués, transformant une légère différence initiale en un écart considérable.

Définition des bornes

Dans de nombreux problèmes, seuls les premiers termes d'un polynôme nous intéressent. Il est alors pratique d'imposer une borne supérieure au degré, c'est-à-dire de travailler modulo \(x^n\). Cette approche possède les propriétés suivantes, garantissant la gestion efficace des séries infinies :

Notation : \([x^n]F(x)\) et \(F[n]\) désignent tous deux le coefficient du terme de degré \(n\) de \(F(x)\).

L'addition et la soustraction de polynômes se définissent de manière évidente :

Inverse modulo pour la multiplication

Deux polynômes \(F(x)\) et \(G(x)\) sont inverses multiplicatifs si \(F(x)*G(x)=1\). Pour calculer \(G(x)\) tel que \(F(x)*G(x) \equiv 1 \pmod{x^n}\), on utilise une approche par dichotomie. Supposons que l'on connaisse un polynôme \(H(x)\) satisfaisant \(F(x)H(x) \equiv 1 \pmod{x^{\lceil n/2 \rceil}}\). L'objectif est de déterminer \(G(x)\) à partir de \(H(x)\).

De la congruence initiale, on déduit que \(G(x)-H(x) \equiv 0 \pmod{x^{\lceil n/2 \rceil}}\). En élevant au carré (en adaptant le module), puis en éliminant le terme \(G(x)^2\) à l'aide de \(F(x)\), on obtient la relation de récurrence :

Voici une implémentation possible, adaptée pour cet article :

#define INITIALIZE_ARRAY(a,n) memset(a,0,sizeof(int)*(n))
#define COPY_ARRAY(dest,src,n) memcpy(dest,src,sizeof(int)*(n))
#define REVERSE_SEGMENT(a,b) reverse(a,b)

int modular_exponentiation(int base, int exp, int modulus) {
    int result = 1;
    base %= modulus;
    while (exp > 0) {
        if (exp & 1) result = (1LL * result * base) % modulus;
        base = (1LL * base * base) % modulus;
        exp >>= 1;
    }
    return result;
}

void modular_reduction(int &value, int modulus) {
    value -= modulus;
    value += (value >> 31) & modulus;
}

int bit_reverse[M], roots[M], transform_limit, log_limit, previous_n = -1;

void initialize_ntt(int n) {
    if (n == previous_n) return;
    previous_n = n;
    for (log_limit = 0, transform_limit = 1; transform_limit <= n; log_limit++, transform_limit <<= 1);
    for (int i = 0; i < transform_limit; i++)
        bit_reverse[i] = (bit_reverse[i >> 1] >> 1) | ((i & 1) << (log_limit - 1));
    for (int len = 1; len < transform_limit; len <<= 1) {
        int w = modular_exponentiation(3, (MOD - 1) / (len << 1), MOD);
        roots[len] = 1;
        for (int j = 1; j < len; j++)
            roots[len + j] = (1LL * roots[len + j - 1] * w) % MOD;
    }
}

void number_theoretic_transform(int *coefficients, int inverse) {
    for (int i = 0; i < transform_limit; i++)
        if (i > bit_reverse[i])
            swap(coefficients[i], coefficients[bit_reverse[i]]);
    for (int len = 1; len < transform_limit; len <<= 1)
        for (int start = 0; start < transform_limit; start += (len << 1))
            for (int k = 0; k < len; k++) {
                int u = coefficients[start + k];
                int v = (1LL * roots[len + k] * coefficients[start + len + k]) % MOD;
                modular_reduction(coefficients[start + k] = u + v, MOD);
                modular_reduction(coefficients[start + len + k] = u + MOD - v, MOD);
            }
    if (inverse) return;
    int inv_n = modular_exponentiation(transform_limit, MOD - 2, MOD);
    for (int i = 0; i < transform_limit; i++)
        coefficients[i] = (1LL * coefficients[i] * inv_n) % MOD;
}

void polynomial_multiply(int n, int m, int *poly1, int *poly2, int *result) {
    initialize_ntt(n + m);
    static int temp1[M], temp2[M];
    COPY_ARRAY(temp1, poly1, n); INITIALIZE_ARRAY(temp1 + n, transform_limit - n);
    COPY_ARRAY(temp2, poly2, m); INITIALIZE_ARRAY(temp2 + m, transform_limit - m);
    number_theoretic_transform(temp1, 1);
    number_theoretic_transform(temp2, 1);
    for (int i = 0; i < transform_limit; i++)
        result[i] = (1LL * temp1[i] * temp2[i]) % MOD;
    REVERSE_SEGMENT(result + 1, result + transform_limit);
    number_theoretic_transform(result, 0);
    INITIALIZE_ARRAY(result + n, transform_limit - n);
}

void polynomial_inverse(int n, int *input, int *inverse) {
    static int temp[M];
    if (n == 1) {
        inverse[0] = modular_exponentiation(input[0], MOD - 2, MOD);
        return;
    }
    polynomial_inverse((n + 1) >> 1, input, inverse);
    initialize_ntt(n << 1);
    COPY_ARRAY(temp, input, n); INITIALIZE_ARRAY(temp + n, transform_limit - n);
    INITIALIZE_ARRAY(inverse + n, transform_limit - n);
    number_theoretic_transform(temp, 1);
    number_theoretic_transform(inverse, 1);
    for (int i = 0; i < transform_limit; i++)
        inverse[i] = (1LL * inverse[i] * (2 + MOD - 1LL * inverse[i] * temp[i] % MOD)) % MOD;
    REVERSE_SEGMENT(inverse + 1, inverse + transform_limit);
    number_theoretic_transform(inverse, 0);
    INITIALIZE_ARRAY(inverse + n, transform_limit - n);
}

Division euclidienne de polynômes

Étant donné \(F(x)\) et \(G(x)\), on cherche le quotient \(Q(x)\) et le reste \(R(x)\) tels que \(F(x)=G(x)Q(x)+R(x)\), où \(\deg(Q)=n-m\) et \(\deg(R)[F^R(x)\equiv G^R(x)Q^R(x)\pmod{x^{n-m+1}}]Ainsi, \(Q^R(x)\) peut être calculé via l'inverse modulo de \(G^R(x)\). Le reste est ensuite obtenu par \(R(x)=F(x)-G(x)Q(x)\).

void polynomial_division(int n, int m, int *dividend, int *divisor, int *quotient, int *remainder) {
    static int rev_dividend[M], rev_divisor[M], temp_inv[M];
    COPY_ARRAY(rev_dividend, dividend, n);
    COPY_ARRAY(rev_divisor, divisor, m);
    REVERSE_SEGMENT(rev_dividend, rev_dividend + n);
    REVERSE_SEGMENT(rev_divisor, rev_divisor + m);
    for (int i = n - m + 1; i < (n << 1); i++) rev_divisor[i] = 0;
    polynomial_inverse(n - m + 1, rev_divisor, temp_inv);
    polynomial_multiply(n, n - m + 1, rev_dividend, temp_inv, quotient);
    for (int i = n - m + 1; i < (n << 1); i++) quotient[i] = 0;
    REVERSE_SEGMENT(quotient, quotient + n - m + 1);
    polynomial_multiply(n - m + 1, m, quotient, divisor, remainder);
}

Interpolation polynomiale de Lagrange -------------------------------------

Étant donné \(n\) couples de points \((x_i, y_i)\), le polynôme d'interpolation \(f(x)\) passant par ces points est :

Forme barycentrique de l'interpolation de Lagrange

En posant \(g=\prod_{i=0}^n (k-x_i)\), la formule se réécrit sous la forme barycentrique :

Une implémentation naïve en \(\mathcal{O}(n^2)\) peut être suffisante pour de petites valeurs de \(n\), comme dans l'exemple ci-dessous :

const int MAX_POINTS = 2e3 + 10;
int n, eval_point, xs[MAX_POINTS], ys[MAX_POINTS];

int main() {
    // Lecture des données...
    long long final_result = 0;
    for (int i = 1; i <= n; i++) {
        long long numerator = 1, denominator = 1;
        for (int j = 1; j <= n; j++) {
            if (i == j) continue;
            numerator = numerator * (eval_point - xs[j] + MOD) % MOD;
            denominator = denominator * (xs[i] - xs[j] + MOD) % MOD;
        }
        final_result = (final_result + numerator * modular_exponentiation(denominator, MOD - 2, MOD) % MOD * ys[i] % MOD) % MOD;
    }
    // Affichage du résultat...
}

Développement de Taylor et itération de Newton ----------------------------------------------

Le développement de Taylor d'une fonction \(F(x)\) autour d'un point \(x_0\) donne :

Logarithme polynômial ---------------------

Le loagrithme formel \(B(x)=\ln(A(x))\) satisfait la relation \(B'(x)=\frac{A'(x)}{A(x)}\). On peut donc le calculer en dérivant \(A(x)\), en multipliant par l'inverse de \(A(x)\), puis en intégrant.

void polynomial_derivative(int n, int *input, int *output) {
    for (int i = 0; i < n - 1; i++) output[i] = (1LL * (i + 1) * input[i + 1]) % MOD;
    output[n - 1] = 0;
}

void polynomial_integral(int n, int *input, int *output) {
    output[0] = 0;
    for (int i = 1; i <= n; i++) output[i] = (1LL * input[i - 1] * inverse_table[i]) % MOD;
}

void polynomial_logarithm(int n, int *input, int *output) {
    static int derivative[M], inverse_poly[M];
    polynomial_derivative(n, input, derivative);
    INITIALIZE_ARRAY(inverse_poly, n);
    polynomial_inverse(n, input, inverse_poly);
    polynomial_multiply(n, n, derivative, inverse_poly, derivative);
    polynomial_integral(n, derivative, output);
}

Exponentielle polynômiale -------------------------

Étant donné un polynôme \(A(x)\) avec \(A(0)=0\), on cherche \(F(x)=\exp(A(x))\). La relation \(G(F(x))=\ln(F(x))-A(x)=0\) permet d'appliquer l'itération de Newton. En utilisant \(G'(F(x))=1/F(x)\), on obtient :

void polynomial_exponential(int n, int *input, int *output) {
    if (n == 1) {
        output[0] = 1;
        return;
    }
    static int temp_log[M];
    polynomial_exponential((n + 1) >> 1, input, output);
    INITIALIZE_ARRAY(temp_log, n);
    polynomial_logarithm(n, output, temp_log);
    for (int i = 0; i < n; i++)
        temp_log[i] = (input[i] + MOD - temp_log[i]) % MOD;
    temp_log[0] = (temp_log[0] + 1) % MOD; // Correction pour l'incrémentation implicite
    polynomial_multiply(n, n, output, temp_log, output);
    INITIALIZE_ARRAY(output + n, transform_limit - n);
}

Racine carrée polynômiale -------------------------

On cherche \(F(x)\) tel que \(F^2(x) \equiv A(x) \pmod{x^n}\), avec \(a_0=1\). En posant \(G(F(x))=F^2(x)-A(x)\), l'itération de Newton donne :

void polynomial_sqrt(int n, int *input, int *output) {
    if (n == 1) {
        output[0] = 1; // Par hypothèse, la racine carrée du premier coefficient est 1
        return;
    }
    static int inverse_sqrt[M], temp[M];
    polynomial_sqrt((n + 1) >> 1, input, output);
    INITIALIZE_ARRAY(inverse_sqrt, n);
    polynomial_inverse(n, output, inverse_sqrt);
    polynomial_multiply(n, n, inverse_sqrt, input, temp);
    int half = modular_exponentiation(2, MOD - 2, MOD);
    for (int i = 0; i < n; i++)
        output[i] = (1LL * (output[i] + temp[i]) * half) % MOD;
}

Exponentiation rapide de polynômes ----------------------------------

Pour calucler \(B(x) \equiv A^k(x) \pmod{x^n}\), on utilise la propriété \(\ln(B(x)) = k \cdot \ln(A(x))\), suivi d'une exponentiation. Pour des exposants très grands, des propriétés modulo liées aux nombres premiers (comme \(a_0=1\)) doivent être prises en compte.

void polynomial_power(int n, int exponent, int *input, int *output) {
    static int log_poly[M];
    polynomial_logarithm(n, input, log_poly);
    for (int i = 0; i < n; i++)
        log_poly[i] = (1LL * log_poly[i] * exponent) % MOD;
    polynomial_exponential(n, log_poly, output);
}

Techniques et optimisations ---------------------------

Convolution avec décalage

Pour calculer \(C_k=\sum_{i=k}^n A_i B_{i-k}\), on peut renverser le second polynôme et effectuer une convolution standard.

Translation d'un polynôme

Pour obtenir les coefficients de \(F(x+c)\) à partir de \(F(x)\), on utilise une transformation factorielle et une convolution.

Conseils d'implémentation

- Pour les petites tailles, une multiplication par force brute peut surpasser la NTT. - L'utilisation de fonctions inline et d'optimisations mémoire (comme memset) accélère le code. - L'initialisation statique des tableaux dans les fonctions récursives évite des allocations coûteuses.

Exercices et applications -------------------------

Calcul de forces électrostatiques

Le calcul d'une force totale impliquant des termes en \(1/r^2\) peut être formulé comme une convolution après inversion des indices.

Optimisation d'un cadeau

Le problème de minimiser la somme des carrés des différences après un ajustement constant peut se décomposer en une optimisation directe et une convolution pour maximiser un produit scalaire, exploitant la structure circulaire.

Sommations et différences itérées

Les opérateurs de somme partielle et de différence itérée correspondent à des convolutions avec des polynômes de la forme \((1-x)^{\pm k}\), calculables via l'exponentiation rapide.

Le flair du pionnier

Ce problème illustre l'importance de la méthode de contribution. Pour chaque élément de la séquence, le nombre de séquences d'intervalles valides peut être exprimé par des coefficients binomiaux, menant à une solution par convolution.

Ressources complémentaires --------------------------

Pour approfondir, consulter des notes de cours sur les polynômes ou des tutoriels spécialisés sur l'algèbre linéaire et les transformations rapides de Fourier.