多项式求逆

前置技能

FFT、NTT、多项式乘法。

原理推理

现在我们有一个这样的式子:

\[ A(x)B(x) \equiv C(x) \ (mod \ x^n) \]

现在我们已知\(B(x), C(x)\)的信息,而我们计算答案的时候需要用\(A(x)\)进行运算。这个时候,我们需要求出\(B(x)\)的逆元。多项式求逆元,我们一般是用倍增的方式来求解。假如我们有\(B(x)\)在模\(x^{\lceil \frac{x}{2} \rceil}\)意义下的逆元\(D(x)\):

\[ B(x)D(x) \equiv 1 \ (mod \ x^{\lceil \frac{x}{2} \rceil}) \]

现在我们需要求\( x^n \)意义下的逆元,我们现在有两个式子:

\[ \begin{cases} B(x)D(x) \equiv 1 \ (mod \ x^{\lceil \frac{x}{2} \rceil}) \\ B(x)B^{-1}(x) \equiv 1 \ (mod \ x^{\lceil \frac{x}{2} \rceil}) \end{cases} \]

上下相减:

\[ B(x)(D(x) – B^{-1}(x)) \equiv 0 \ (mod \ x^{\lceil \frac{x}{2} \rceil}) \\ D(x) – B^{-1}(x) \equiv 0 \  (mod \ x^{\lceil \frac{x}{2} \rceil}) \]

同时平方:

\[ D^2(x) + B{-2}(x) – 2D(x)B^{-1}(x) \equiv 0 \ (mod \ n) \]

其中,这里为什么模意义突然变成了\(n\)的原因如下:

我可能会这样讲:根据\(B(x) – B^{-1}(x) \equiv 0 \pmod {x^{\lceil n/2 \rceil}}\),左边一定是一个形如\(k\cdot x^{\lceil n/2 \rceil}\)的多项式,它的平方(即\((B(x) – B^{-1}(x))^2\))为\(k^2 x^n\)(\(n\)为偶数)或\(k^2 x^{n+1}\)(\(n\)为奇数) ,显然二者均能被\(x_n\)整除

src: http://blog.miskcoo.com/2015/05/polynomial-inverse, PLANET6174

剩下的这个乘法用 NTT 解决就行了。

代码

// P3803.cpp
// NTT version;
#include <bits/stdc++.h>
#define ll long long

using namespace std;

const int mod = 998244353, G = 3, Gi = 332748118, MAX_N = 3e6 + 2000;

int n, m, rev[MAX_N], mx_bit, mx_pow;

ll ai[MAX_N], bi[MAX_N];

ll quick_pow(ll bas, ll tim)
{
    ll ans = 1;
    bas %= mod;
    while (tim)
    {
        if (tim & 1)
            ans = ans * bas % mod;
        bas = bas * bas % mod;
        tim >>= 1;
    }
    return ans;
}

void ntt_initialize()
{
    mx_pow = 2, mx_bit = 1;
    while (n + m >= mx_pow)
        mx_pow <<= 1, mx_bit++;
    for (int i = 0; i < mx_pow; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (mx_bit - 1));
}

inline void ntt(ll *arr, int dft)
{
    for (int i = 0; i < mx_pow; i++)
        if (i < rev[i])
            swap(arr[i], arr[rev[i]]);
    for (int step = 1; step < mx_pow; step <<= 1)
    {
        ll omega_n = quick_pow(dft == 1 ? G : Gi, (mod - 1) / (step << 1));
        for (int j = 0; j < mx_pow; j += (step << 1))
        {
            ll omega_nk = 1;
            for (int k = j; k < j + step; k++, omega_nk = (omega_nk * omega_n) % mod)
            {
                int x = arr[k], y = omega_nk * arr[k + step] % mod;
                arr[k] = (x + y) % mod;
                arr[k + step] = (x - y + mod) % mod;
            }
        }
    }
}

int main()
{
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; i++)
        scanf("%lld", &ai[i]), ai[i] = (ai[i] + mod) % mod;
    for (int i = 0; i <= m; i++)
        scanf("%lld", &bi[i]), bi[i] = (bi[i] + mod) % mod;
    ntt_initialize();
    ntt(ai, 1), ntt(bi, 1);
    for (int i = 0; i < mx_pow; i++)
        ai[i] = (ai[i] * bi[i] % mod);
    ntt(ai, -1);
    ll inv = quick_pow(mx_pow, mod - 2);
    for (int i = 0; i <= n + m; i++)
        printf("%d ", ai[i] * inv % mod);
    return 0;
}

发布者

Kal0rona

江西师大附中全机房最弱

发表评论

电子邮件地址不会被公开。 必填项已用*标注