@总结 - 3@ 多项式的进阶运算 —— 多项式求逆，除法与取余

最新推荐文章于 2021-04-27 20:56:17 发布

weixin_30949361

最新推荐文章于 2021-04-27 20:56:17 发布

阅读量346

点赞数

原文链接：http://www.cnblogs.com/Tiw-Air-OAO/p/10187501.html

版权

@0 - 参考资料@
@1 - 前置技能@
@2 - 一些概念@
@3 - 多项式求逆元@
@4 - 多项式除法与取余@

@0 - 参考资料@

Miskcoo's Space 的讲解 1
Miskcoo's Space 的讲解 2

@1 - 前置技能@

这里是本蒟蒻的 FFT 讲解
如果各位 dalao 觉得我这个蒟蒻写得实在不行，可以去看看上面我提到的那位 dalao 的 FFT 讲解 qwq。

@2 - 一些概念@

多项式的度：对于多项式 \(A(x)\)，我们称它的最高次数为 \(degA\)。

多项式的除法：对于多项式 \(A(x)\) 与 \(B(x)\)，如果有:
\[A(x) = B(x)*C(x) + R(x)\]
则称 \(C(x)\) 是 \(A(x)\) 除以 \(B(x)\) 的商，\(R(x)\) 是 \(A(x)\) 除以 \(B(x)\) 的余数。

类比整数，我们可以得到多项式的 “同余” ：
\[A(x)\equiv R(x)\mod B(x)\]

多项式的逆元：对于多项式 \(A(x)\)，如果存在多项式 \(B(x)\)，使得：
\[A(x)*B(x)\equiv 1\mod x^{n}\]

则称 \(B(x)\) 为 \(A(x)\) 模 \(x^{n}\) 意义下的逆元，记作 \(A^{-1}(x)\)。

@3 - 多项式求逆元@

@理论推导@

考虑用倍增法来解决这一问题。

当 \(n = 1\) 时，\(A(x)\mod x\) 是个常数 \(c\)，此时 \(A^{-1}(x)\equiv c^{-1}\mod x\)。

假如我们已知 \(A^{-1}(x)\equiv B'(x)\mod x^{\lceil\frac{n}{2}\rceil}\)（为什么是向上取整？等会儿你就知道了 qwq）。
怎么求出 \(A^{-1}(x)\equiv B(x)\mod x^n\) 的值呢？

首先可以得到：
\[A(x)*B(x)\equiv 1 \mod x^{\lceil\frac{n}{2}\rceil}\]
两边同时乘 \(B'(x)\)，得到：
\[B'(x)*A(x)*B(x)\equiv B(x)\equiv B'(x) \mod x^{\lceil\frac{n}{2}\rceil}\]
（因为 \(B'(x)\) 是 \(A(x)\) 的逆元）
所以有：
\[B(x)-B'(x) \equiv 0 \mod x^{\lceil\frac{n}{2}\rceil}\]

然后这这一步可能是整个推导里面比较难懂的一部分：
我们把等式两边连同模数平方，得到：
\[(B(x)-B'(x))^2 \equiv 0 \mod x^n\]

为什么它成立呢？
因为 \(B(x)-B'(x) \mod x^{\lceil\frac{n}{2}\rceil} = 0\) ，所以 \(B(x)-B'(x)\) 这个多项式中 \(x^{0\dots \lceil\frac{n}{2}\rceil-1}\) 的系数都等于 0，所以次数最小的系数可能不为 0 的项为 \(x^{0\dots \lceil\frac{n}{2}\rceil}\)。
平方过后，次数最小的系数可能不为 0 的项为 \(x^{\lceil\frac{n}{2}\rceil*2}\)。
然后它就成立了。

这个时候我们就能解释为什么要向上取整而不是向下取整了：向下取整过后可能\(\lfloor\frac{n}{2}\rfloor*2 < n\)，上面的操作就不能够正常完成了。

接上我们刚刚的推导，把括号展开：
\[B(x)^2-2B(x)B'(x)+B'(x)^2 \equiv 0 \mod x^n\]

两边同时乘 \(A(x)\)：
\[B(x)-2B'(x)+B'(x)^2*A(x) \equiv 0 \mod x^n\]

移项：
\[B(x) \equiv 2B'(x)-B'(x)^2*A(x) \mod x^n\]

用 FFT 搞一下多项式乘法。
然后就没了。

总时间复杂度 \(O(n\log n)+O(\frac{n}{2}\log \frac{n}{2}) + O(\frac{n}{4}\log \frac{n}{4}) + ... = O(n\log n)\)。
不要问我为什么时间复杂度是这样，我也不会证明。

@参考代码@

本代码为 luogu4238 的 AC 代码。

#include<cstdio>
#include<algorithm>
using namespace std;
const int G = 3;
const int MOD = 998244353;
const int MAXN = 400000;
int pow_mod(int b, int p) {
    int ret = 1;
    while( p ) {
        if( p & 1 ) ret = 1LL*ret*b%MOD;
        b = 1LL*b*b%MOD;
        p >>= 1;
    }
    return ret;
}
void NTT(int A[], int n, int type) {
    for(int i=0,j=0;i<n;i++) {
        if( i < j ) swap(A[i], A[j]);
        for(int l=(n>>1);(j^=l)<l;l>>=1);
    }
    for(int s=2;s<=n;s<<=1) {
        int t = (s>>1);
        int u = (type == 1) ? pow_mod(G, (MOD-1)/s) : pow_mod(G, (MOD-1)-(MOD-1)/s);
        for(int i=0;i<n;i+=s) {
            for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {
                int x = A[i+j], y = 1LL*A[i+j+t]*p%MOD;
                A[i+j] = (x + y)%MOD;
                A[i+j+t] = (x - y + MOD)%MOD;
            }
        }
    }
    if( type == -1 ) {
        int inv = pow_mod(n, MOD-2);
        for(int i=0;i<n;i++)
            A[i] = 1LL*A[i]*inv%MOD;
    }
}
int tmp[MAXN + 5];
void poly_inv(int *A, int *B, int n) {
    if( n == 1 ) {
        B[0] = pow_mod(A[0], MOD-2);
        return ;
    }
    poly_inv(A, B, (n + 1)>>1);//注意这里必须扩展到 2n 的长度
    int m = 1; while( m < (n<<1) ) m <<= 1;
    for(int i=0;i<n;i++)
        tmp[i] = A[i];
    for(int i=n;i<m;i++)
        tmp[i] = 0;
    NTT(tmp, m, 1); NTT(B, m, 1);
    for(int i=0;i<m;i++)
        B[i] = 1LL*B[i]*(2LL + MOD - 1LL*tmp[i]*B[i]%MOD)%MOD;
    NTT(B, m, -1);
    for(int i=n;i<m;i++)
        B[i] = 0;
}
int A[MAXN + 5], B[MAXN + 5];
int main() {
    int n; scanf("%d", &n);
    for(int i=0;i<n;i++)
        scanf("%d", &A[i]);
    poly_inv(A, B, n);
    for(int i=0;i<n;i++)
        printf("%d%c", B[i], (i == n-1) ? '\n' : ' ');
}

@例题与应用@

@伯努利数@

建议大家可以去看看这一位 dalao 的博客。

定义伯努利数 \(B_i\)，满足：
\[\begin{cases}B_0=1\\ C_n^0B_0+C_n^1B_1+\dots+C_n^nB_n=B_n&n > 1或n = 0\end{cases}\]
可以发现它有递推公式，但是不是很重要，因为用递推公式求伯努利数是 \(O(n^2)\)。

指数型生成函数为 \(\hat F(x) = \dfrac{F_0}{0!} + \dfrac{F_1}{1!}x + \dfrac{F_2}{2!}x^2 + \dots\)

令两个指数型生成函数 \(\hat F(x),\hat G(x)\) 相乘得到新生成函数 \(\hat H(x)\)，有：
\[H_i=\sum_{j+k=i}C_{j+k}^{k}*F(j)*G(k)\]

这个你只需要暴力展开验证即可。

再来看伯努利数的性质的第二条，等式左边很像是我们生成函数的乘法。事实上，我们的确可以把左边看成伯努利数的指数型生成函数与系数全为 1 的数列（\(e^x\)）的指数型生成函数作乘法。

也就是说：\[\hat B(x)e^x=\hat B(x) + x\]

为什么右边会多出来个 \(+ x\) 呢？注意我们的第二条性质当 \(n = 1\) 时是不成立的。
所以我们得到伯努利数的指数型生成函数： \(\hat B(x) = \dfrac{x}{e^x-1}\)。

将分母泰勒展开，就可以得到： \(\hat B(x) = \dfrac{1}{1+\frac{1}{2!}x+\frac{1}{3!}x^2+\dots}\)。
一个多项式的逆元，可以用我们刚刚介绍的方法来求解了。

那么这玩意儿有什么用呢？它可以用来求自然数的幂和。即求：
\[S(n, k) = 0^k + 1^k+2^k+3^k+\dots+n^k\]
的值（我们钦定 \(0^0 = 1\)）。

将 \(e^{kx}\) 泰勒展开：
\[e^{kx}=1+kx+\dfrac{k^2}{2!}x^2+\dfrac{k^3}{3!}x^3\dots\]

则 \(e^{0x}+e^{1x}+e^{2x}+...+e^{nx}\) 的第 \(k\) 项的值就是 \(S(n, k)\)。

等比数列求和：\(\hat S(x) = \dfrac{e^{(n+1)x}-1}{e^x-1}\)（注意左边是 \(S\) 的生成函数）。
注意到分母没有？我们可以把伯努利数弄过来了。
\[\hat S(x) = \hat B(x)\dfrac{e^{(n+1)x}-1}{x}=\hat B(x)(n+1+\dfrac{(n+1)^2}{2!}x+\dots)\]

又是两个指数型生成函数乘积，其中另一个数列为 \(\{\dfrac{n+1}{1},\dfrac{(n+1)^2}{2},\dfrac{(n+1)^3}{3}\dots\}\)。

所以：
\[S_i=\sum_{j+k=i}C_{j+k}^{k}B_j\dfrac{(n+1)^{k+1}}{k+1}\]

啥？你说这个公式和网上的其他公式不一样？你把组合数拆成阶乘，你再来看看它们是不是等价的？

总结一下，关于伯努利数的推导，你需要反复用生成函数，泰勒展开，生成函数，泰勒展开……
然后你只需要记住 \(\hat B(x) = \dfrac{1}{\sum_{i=0}\frac{1}{(i+1)!}x^i}\) 与 \(S(n,i)=\sum_{j+k=i}C_{j+k}^{k}B_j\dfrac{(n+1)^{k+1}}{k+1}\)。

@4 - 多项式除法与取余@

@理论推导@

\[A(x) = B(x)*C(x) + R(x)\]
其中 \(degA = n\)， \(degB = m\)，\(degC = n-m\)，\(degR < m\)。
我们已知 \(A(x)\) 与 \(B(x)\)，需要求解 \(C(x)\) 与 \(R(x)\)。

如果没有 \(R(x)\) 多好啊，我们就可以直接在两边乘上 \(B(x)\) 的逆元了。
考虑怎么消除 \(R(x)\) 的影响。

我们假设有一个多项式 \(P(x) = p_0 + p_1x + p_2x^2+...+p_kx^{k}\)
将 \(x^{-1}\) 代入该多项式得到：
\[P(x^{-1}) = p_0 + p_1x^{-1}+p_2x^{-2}+...+p_kx^{-k}\]
再乘 \(x^k\) 得到：
\[x^kP(x^{-1}) = p_0x^k + p_1x^{k-1}+p_2x^{k-2}+...+p_k\\ =p_k+p_{k-1}x+...+p_0x^k\]

对于初始的 P 来说相当于把系数翻转，因此记 \(P^r(x) = x^kP(x^{-1})\)。
【这种操作……我真的想不到 TAT】

好，我们回到初始的除法式子。
\[A(x) = B(x)*C(x) + R(x)\]
代入 \(x^{-1}\) ，再两边同时乘 \(x^n\)：
\[x^nA(x^{-1}) = x^mB(x^{-1})*x^{n-m}C(x^{-1}) + x^{m-1}R(x^{-1})*x^{n-m+1}\]
\[A^r(x)=B^r(x)*C^r(x)+x^{n-m+1}R^r(x)\]
这样余式 \(R(x)\) 的最低次数就变成 \(x^{n-m+1}\) 啦。

两边同时取模，得到：
\[A^r(x) \equiv B^r(x)*C^r(x)\mod x^{n-m+1}\]
然后就会发现我们达到了我们最初的目的：消去 \(R(x)\) 的影响。
两边乘一下逆元 \((B^r(x))^{-1}\)就 OK。

因为多项式翻转过后，它的 \(deg\) 不会变化，所以 \(degC = n-m\)。
所以我们可以直接解出来 \(C^r(x)\) 的值。

最后余式 \(R(x)\) 直接用 \(A(x) - B(x)*C(x)\) 即可~

@参考代码@

本代码为 luogu4512 的 AC 代码。

#include<cstdio>
#include<algorithm>
using namespace std;
const int G = 3;
const int MOD = 998244353;
const int MAXN = 400000;
int pow_mod(int b, int p) {
    int ret = 1;
    while( p ) {
        if( p & 1 ) ret = 1LL*ret*b%MOD;
        b = 1LL*b*b%MOD;
        p >>= 1;
    }
    return ret;
}
struct Polynomial{
    int tmp[MAXN + 5], tmp1[MAXN + 5], tmp2[MAXN + 5], tmp3[MAXN + 5];
    void poly_copy(int *A, int *B, int n) {
        for(int i=0;i<n;i++)
            A[i] = B[i];
    }
    void poly_clear(int *A, int l, int r) {
        for(int i=l;i<r;i++)
            A[i] = 0;
    }
    void poly_revcopy(int *A, int *B, int n) {
        poly_copy(A, B, n);
        for(int i=0,j=n-1;i<j;i++,j--)
            swap(A[i], A[j]);
    }
    void ntt(int *A, int n, int type) {
        for(int i=0, j=0;i<n;i++) {
            if( i < j ) swap(A[i], A[j]);
            for(int l=(n>>1);(j^=l)<l;l>>=1);
        }
        for(int s=2;s<=n;s<<=1) {
            int t = (s>>1);
            int u = (type == -1) ? pow_mod(G, (MOD-1) - (MOD-1)/s) : pow_mod(G, (MOD-1)/s);
            for(int i=0;i<n;i+=s) {
                for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {
                    int x = A[i+j], y = 1LL*p*A[i+j+t]%MOD;
                    A[i+j] = (x + y)%MOD, A[i+j+t] = (x + MOD - y)%MOD;
                }
            }
        }
        if( type == -1 ) {
            int inv = pow_mod(n, MOD-2);
            for(int i=0;i<n;i++)
                A[i] = 1LL*A[i]*inv%MOD;
        }
    }
    void poly_mul(int *A, int *B, int *C, int n, int m) {
        int len; for(len = 1;len < n + m - 1;len <<= 1);
        ntt(A, len, 1); ntt(B, len, 1);
        for(int i=0;i<len;i++)
            C[i] = 1LL*A[i]*B[i]%MOD;
        ntt(C, len, -1);
    }
    void poly_inv(int *A, int *B, int n) {
        if( n == 1 ) {
            B[0] = pow_mod(A[0], MOD-2);
            return ;
        }
        poly_inv(A, B, (n + 1) >> 1);
        int len; for(len = 1;len < (n<<1);len <<= 1);//这里必须要扩展到 2n
        poly_copy(tmp, A, n); poly_clear(tmp, n, len);
        ntt(tmp, len, 1); ntt(B, len, 1);
        for(int i=0;i<len;i++)
            B[i] = 1LL*B[i]*(2 + MOD - 1LL*tmp[i]*B[i]%MOD)%MOD;
        ntt(B, len, -1); poly_clear(B, n, len);
    }
    //我有着良好的压行习惯 OuO 
    void poly_divide(int *A, int *B, int *C, int *R, int n, int m) {
        int len; for(len = 1;len < (n<<1);len <<= 1);//这里必须要扩展到 2n
        
        poly_revcopy(tmp1, B, m); poly_inv(tmp1, tmp2, n-m+1);
        poly_revcopy(tmp1, A, n);
        ntt(tmp1, len, 1); ntt(tmp2, len, 1);
        for(int i=0;i<len;i++)
            tmp3[i] = 1LL*tmp1[i]*tmp2[i]%MOD;
        ntt(tmp3, len, -1); poly_revcopy(C, tmp3, n-m+1);
            
        poly_copy(tmp1, B, len); ntt(tmp1, len, 1);
        poly_copy(tmp2, C, len); ntt(tmp2, len, 1);
        for(int i=0;i<len;i++) tmp3[i] = 1LL*tmp1[i]*tmp2[i]%MOD;
        ntt(tmp3, len, -1);
        for(int i=0;i<len;i++) R[i] = (A[i] + MOD - tmp3[i])%MOD;
    }
}poly;
int f[MAXN + 5], g[MAXN + 5];
int r[MAXN + 5], q[MAXN + 5];
int main() {
    int n, m;
    scanf("%d%d", &n, &m); n++, m++;
    for(int i=0;i<n;i++)
        scanf("%d", &f[i]);
    for(int i=0;i<m;i++)
        scanf("%d", &g[i]);
    poly.poly_divide(f, g, q, r, n, m);
    for(int i=0;i<n-m+1;i++)
        printf("%d ", q[i]);
    puts("");
    for(int i=0;i<m-1;i++)
        printf("%d ", r[i]);
}