526互联

多项式全家桶

发布时间 2023-07-10 18:42:50作者: crimson000

多项式基本运算

这个博客主要用来放一些多项式的运算的板子，大部分都来自于洛谷。

多项式乘法

NTT 求个卷积即可。


#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N = 3e6 + 10, mod = 998244353, GG = 3, Gi = 332748118;
ll qmi(ll a, ll k, ll p)
{
    ll res = 1;
    while (k)
    {
        if (k & 1)
            res = res * a % p;
        a = a * a % p;
        k >>= 1;
    }
    return res;
}
int F[N], G[N];
int rev[N];
int n, m;
int lim = 1, len;

void NTT(int a[], int opt)
{
    for (int i = 0; i < lim; i++)
        if (i < rev[i])
            swap(a[i], a[rev[i]]);
    int up = log2(lim);
    for (int dep = 1; dep <= up; dep++)
    {
        int m = 1 << dep;
        int gn;
        if(opt == 1) gn = qmi(GG, (mod - 1) / m, mod);
        else gn = qmi(Gi, (mod - 1) / m, mod);
        for (int k = 0; k < lim; k += m)
        {
            int g = 1;
            for (int j = 0; j < m / 2; j ++)
            {
                int t = (ll)g * a[j + k + m / 2] % mod;
                int u = a[j + k];
                a[j + k] = (u + t) % mod;
                a[j + k + m / 2] = (u - t + mod) % mod;
                g = (ll)g * gn % mod;
            }
        }
    }
    if (opt == -1)
    {
        int inv = qmi(lim, mod - 2, mod);
        for (int i = 0; i < lim; i ++)
            a[i] = (ll)a[i] * inv % mod;
    }
}

int main()
{
    n = read(), m = read();
    for (int i = 0; i <= n; i ++)
        F[i] = (read() + mod) % mod;
    for (int i = 0; i <= m; i ++)
        G[i] = (read() + mod) % mod;

    while (lim <= n + m)
        lim <<= 1, len ++;
    for (int i = 0; i < lim; i ++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));

    NTT(F, 1), NTT(G, 1);
    for (int i = 0; i <= lim; i ++)
        F[i] = (ll)F[i] * G[i] % mod;
    NTT(F, -1);

    for (int i = 0; i <= n + m; i ++)
        cout << F[i] % mod << ' ';

    return 0;
}

多项式求逆

原理

\[\begin{aligned}A\times B &\equiv 1\ (mod \ x^n)\\ A\times B' &\equiv 1\ (mod\ x^{\frac{n}{2}}) \\A \times (B-B')&\equiv 0 \ (mod \ x^{\frac{n}{2} }) \\(B-B')^2&\equiv 0\ (mod \ x^n) \\B^2-2BB'+B'^2&\equiv 0 \ (mod\ x^n)\\A(B^2-2BB'+B'^2 )&\equiv 0 \ (mod\ x^n) \\B-2B'+AB'^2&\equiv 0\ (mod\ x^n)\\ B&\equiv 2B'-AB'^2\ (mod\ x^n)\end{aligned} \]

代码：

（警示后人：慎用memset）

inline void mul(int a[], int b[], int to[], int lim)
{
    for(int i = (lim >> 1); i <= lim; i ++ ) X[i] = Y[i] = 0;
    for(int i = 0; i < (lim >> 1); i ++ )
        X[i] = a[i] % P, Y[i] = b[i] % P;
    NTT(X, lim, 1), NTT(Y, lim, 1);
    for(int i = 0; i < lim; i ++ ) 
        to[i] = (ll)X[i] * Y[i] % P;
    NTT(to, lim, -1);
}

ll b[2][N];
void inv(ll a[], ll to[], ll n)
{
    int cur = 0;
    b[cur][0] = qmi(a[0], P - 2, P);
    int base = 1, lim = 2, len = 1;
    calc(lim, len);
    while(base <= (n + n))
    {
        cur ^= 1;
        memset(b[cur], 0, sizeof b[cur]);
        for(int i = 0; i < base; i ++ ) b[cur][i] = b[cur ^ 1][i] * 2 % P;
        mul(b[cur ^ 1], b[cur ^ 1], b[cur ^ 1], lim);
        mul(b[cur ^ 1], a, b[cur ^ 1], lim);
        for(int i = 0; i < base; i ++ )
            b[cur][i] = (b[cur][i] - b[cur ^ 1][i] + P) % P;
        base <<= 1, lim <<= 1, len ++;
        calc(lim, len);
    }
    for(int i = 0; i < lim; i ++ ) to[i] = b[cur][i];
}

MTT（任意模数多项式乘法）

取 3 个模数跑 9 遍 NTT，最后中国剩余定理合并。

#define LOCAL
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N = 4e5 + 10, G1 = 998244353, G2 = 1004535809, G3 = 469762049, G = 3;
int n, m, p;
int lim = 1, len;
int rev[N];
int A[3][N], B[3][N];
int a1[N], a2[N], a3[N];
ll ans[N];

__int128 qmi(__int128 a, __int128 k, __int128 p)
{
    __int128 res = 1;
    while(k)
    {
        if(k & 1) res = res * a % p;
        a = a * a % p;
        k >>= 1;
    }
    return res;
}

__int128 inv(__int128 a, __int128 mod)
{
    return qmi(a, mod - 2, mod);
}

void NTT(int a[], int mod, int opt)
{
    for(int i = 0; i < lim; i ++ )
        if(i < rev[i])
            swap(a[i], a[rev[i]]);
    int up = log2(lim);
    for(int dep = 1; dep <= up; dep ++ )
    {
        int m = 1 << dep;
        int gn;
        if(opt == 1) gn = qmi(G, (mod - 1) / m, mod);
        else gn = qmi(qmi(G, mod - 2, mod), (mod - 1) / m, mod);
        for(int j = 0; j < lim; j += m)
        {
            int g = 1;
            for(int k = 0; k < m / 2; k ++ )
            {
                int t = (ll)a[j + k + m / 2] * g % mod;
                int u = a[j + k];
                a[j + k] = (u + t) % mod;
                a[j + k + m / 2] = ((ll)u - t + mod) % mod;
                g = (ll)gn * g % mod;
            }
        }
    }
    if(opt == -1)
    {
        int inv = qmi(lim, mod - 2, mod);
        for(int i = 0; i < lim; i ++ ) a[i] = (ll)a[i] * inv % mod;
    }
}

int main()
{
    n = read(), m = read(), p = read();
    for(int i = 0; i <= n; i ++ ) A[0][i] = read(), A[1][i] = A[2][i] = A[0][i];
    for(int i = 0; i <= m; i ++ ) B[0][i] = read(), B[1][i] = B[2][i] = B[0][i];

    while(lim <= n + m + 2) lim <<= 1, len ++;
    for(int i = 0; i < lim; i ++ )
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));

    NTT(A[0], G1, 1), NTT(B[0], G1, 1);
    for(int i = 0; i < lim; i ++ ) a1[i] = (ll)A[0][i] * B[0][i] % G1;
    NTT(a1, G1, -1);

    NTT(A[1], G2, 1), NTT(B[1], G2, 1);
    for(int i = 0; i < lim; i ++ ) a2[i] = (ll)A[1][i] * B[1][i] % G2;
    NTT(a2, G2, -1);

    NTT(A[2], G3, 1), NTT(B[2], G3, 1);
    for(int i = 0; i < lim; i ++ ) a3[i] = (ll)A[2][i] * B[2][i] % G3;
    NTT(a3, G3, -1);

    for(int i = 0; i < lim; i ++ )
    {
        __int128 M = (__int128)G1 * G2 * G3;
        ans[i] = ((__int128)a1[i] * M / G1 * inv(M / G1, G1) % M
                    + (__int128)a2[i] * M / G2 * inv(M / G2, G2) % M + 
                    (__int128)a3[i] * M / G3 * inv(M / G3, G3) % M) % M % p;
    }

    for(int i = 0; i <= n + m; i ++ ) cout << ans[i] << ' ';
    
    return 0;
}

多项式求导/积分

求导公式：

\[(\sum \limits _{i=0}^na_ix^i)'=\sum \limits _{i=0}^{n-1}a_{i+1}(i+1)x^i \]

void derivative(int a[])
{
    for(int i = 0; i < n - 1; i ++ )
        d[i] = (ll)(i + 1) * a[i + 1] % P;
    d[n - 1] = 0;
}

积分公式：

\[\int (\sum \limits _{i=0}^na_ix^i)=\sum \limits _{i=1}^{n+1}\frac{a_{i-1}}{i} x^i \]

void integral(int a[])
{
    for(int i = 1; i < n; i ++ )
        in[i] = (ll)qmi(i, P - 2, P) * a[i - 1] % P;
    in[0] = 0;
}

多项式对数函数

推导：

\[\begin{aligned} g(x)&=\ln (A(x)) \\ &=f(A(x)) \end{aligned} \]

两边同时求导，根据求导公式 \(f(g(x))'=f'(g(x))g'(x)\) 可得

\[\begin{aligned} g'(x) &=f(A(x)) \\ &=\ln (A(x))'A'(x)\\ &=\frac{A'(x)}{A(x)} \end{aligned} \]

两边再同时积分，即可得 \(g(x)=\int \frac{A'(x)}{A(x)}\)。

代码：

ll b[2][N];
void inv(ll a[], ll to[], ll n)
{
    int cur = 0;
    b[cur][0] = qmi(a[0], P - 2, P);
    int base = 1, lim = 2, len = 1;
    calc(lim, len);
    while(base <= (n + n))
    {
        cur ^= 1;
        memset(b[cur], 0, sizeof b[cur]);
        for(int i = 0; i < base; i ++ ) b[cur][i] = b[cur ^ 1][i] * 2 % P;
        mul(b[cur ^ 1], b[cur ^ 1], b[cur ^ 1], lim);
        mul(b[cur ^ 1], a, b[cur ^ 1], lim);
        for(int i = 0; i < base; i ++ )
            b[cur][i] = (b[cur][i] - b[cur ^ 1][i] + P) % P;
        base <<= 1, lim <<= 1, len ++;
        calc(lim, len);
    }
    for(int i = 0; i < lim; i ++ ) to[i] = b[cur][i];
}

void derivative(ll a[], ll to[], ll n)
{
    for(int i = 0; i < n - 1; i ++ )
        to[i] = (ll)(i + 1) * a[i + 1] % P;
    to[n - 1] = 0;
}

void integral(ll a[], ll to[], ll n)
{
    for(int i = 1; i < n; i ++ )
        to[i] = qmi(i, P - 2, P) * a[i - 1] % P;
    to[0] = 0;
}

ll inte[N], INV[N], d[N], LN[N];
void ln(ll a[], ll to[], ll n)
{
    derivative(a, d, n);
    inv(a, INV, n);
    mul(d, INV, d, n, n);
    integral(d, to, n);
}

多项式 exp

牛顿迭代法：\(x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}\)。通过牛顿迭代法，可以得到一些方程的较精确解。

在多项式中，我们也可以运用牛顿迭代法。

\[G(x)=G_0(x)-\frac{F(G_0(x))}{F'(G_0(x))} \]

迭代一次精度就会翻倍，如果 \(F(G_0(x))\equiv 0\ ( mod\ x^{\frac{n}{2} })\)，那么迭代一次会变成 \(F(G(x))\equiv 0\ ( mod\ x^{n })\)。

用牛顿迭代法求解 \(B(x)\equiv e^{A(x)}\)：

\[\begin{aligned} B(x) &\equiv e^{A(x)} \\ \ln B(x) &\equiv A(x) \\ \ln B(x) - A(x)&\equiv 0 \end{aligned} \]

令 \(F(G(x))=\ln G(x) - A(x) =0\)，利用牛顿迭代法：

\[\begin{aligned} G(x)&=G_0(x)-\frac{F(G_0(x))}{F(G_0(x))'} \\ &=G_0(x)-\frac{\ln G_0(x)-A(x)}{\frac{G_0'(x)}{G_0(x)} } \\ &=\frac{G_0(x)(1-\ln G_0(x)+A(x))}{G_0'(x)} \end{aligned} \]

时间复杂度 \(O(n\log n)\)。

PS：在这里多项式求逆我被卡了好久

ll c[N], E[N];
inline void exp(ll a[], ll b[], ll n)
{
    if(n == 1)
    {
        b[0] = 1;
        return;
    }
    exp(a, b, (n + 1) >> 1);
    ln(b, LN, n);
    int lim = 1;
    while(lim <= n + n) lim <<= 1;
    for(int i = 0; i < n; i ++ ) c[i] = ((ll)a[i] - LN[i] + P) % P;
    for(int i = n; i < lim; i ++ ) LN[i] = c[i] = 0;
    c[0] ++;
    mul(c, b, b, n, n);
    for(int i = n; i < lim; i ++ ) b[i] = 0;
}