多项式学习笔记

· · 算法·理论

多项式乘法本质上是加法卷积

[F(x)*G(x)]_i=x^i\sum_{k=0}^{i}F(i)G(k-i)

FFT

$\mathrm {FFT}$:分为$\mathrm {DFT}$和$\mathrm {IDFT}

下面为了方便,设n2的正整数次幂

$$ A_1(x)=a_1x+a_3x^3+a_5x^5+\dots+a_{n-1}x^{n-1} \\ A_2(x)=a_0+a_2x^2+a_4x^4+\dots+a_{n}x^{n} $$ 有 $$ k \in [0,\frac{n}{2}] , A(\omega_{n}^{k})=A_1(\omega_{\frac{n}{2}}^{k})+ \omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^{k}) \\ k\in (\frac{n}{2},n] , A(\omega_{n}^{k})=A_1(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^{k}) \\ $$ $\mathrm{IDFT}$:用$\mathrm {DFT}$中求得的$n$个点值$(\omega_{n}^{k},A(\omega_{n}^{k}))$还原多项式的系数,设$x^k$的系数为$c_k$,记$y_k=A(\omega_{n}^{k})$,有 $$ c_k=\frac{1}{n}\sum_{i=0}^{n-1}y_i(\omega_{n}^{-k})^i $$ 设多项式$B(x)=y_0+y_1x^1+y_2x^2+\dots+y_{n-1}x^{n-1}$,那么$c_k=\frac{1}{n}B(\omega_{n}^{-k})

再做一次\mathrm{DFT}即可

实现:

迭代求解

发现按照偶数项放在左边,奇数项放在右边求解时,初始数组中二进制表示为s的下标的点到达最下面一层时其二进制表示正好是s翻转后的二进制表示,预处理rev_s表示s翻转后的值,做\mathrm {FFT}时将rev_s翻上来。

枚举mid=2^k表示当前处理长度为2\times mid的区间,将这个区间分为左边mid个和右边mid个归并答案即可

基础代码:

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const int N=300010;
const double PI=acos(-1);
int n,m;
struct Complex { double x,y; };
Complex operator + (Complex a,Complex b) { return {a.x+b.x,a.y+b.y}; }
Complex operator - (Complex a,Complex b) { return {a.x-b.x,a.y-b.y}; }
Complex operator * (Complex a,Complex b) { return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; }
Complex a[N],b[N];
int rev[N],bit,tot;

void FFT(Complex a[],int dir)
{
    for (int i=0;i<tot;i++)
    {
        if (i<rev[i])
        {
            swap(a[i],a[rev[i]]);
        }
    }
    for (int mid=1;mid<tot;mid<<=1)
    {
        //处理长度为2mid
        Complex wn1={cos(PI/mid),dir*sin(PI/mid)};
        for (int i=0;i<tot;i+=2*mid)
        {
            Complex wnk={1,0};
            for (int j=0;j<mid;j++,wnk=wnk*wn1)
            {
                Complex x=a[i+j],y=wnk*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;
            }
        }
    }
}

int main()
{
    cin >> n >> m;
    for (int i=0;i<=n;i++) cin >> a[i].x;
    for (int i=0;i<=m;i++) cin >> b[i].x;

    while ((1<<bit)<n+m+1) bit++;
    tot=(1<<bit);

    for (int s=0;s<tot;s++) rev[s]=(rev[s>>1]>>1)|((s&1)<<(bit-1)); //s翻转

    FFT(a,1),FFT(b,1);
    for (int i=0;i<tot;i++) a[i]=a[i]*b[i];
    FFT(a,-1);

    for (int i=0;i<=n+m;i++)
    {
        double c=a[i].x/tot+0.5;
        printf("%d ",(int)c);
    }

    return 0;
}

NTT

在模意义下快速实现多项式乘法的一种方法

你说得对,但是原根是一个数学符号。设 m 是正整数,a 是整数, 若 am 的阶等于 \varphi(m),则称 a 为模 m 的一个原根。假设一个数 gP 的原根,那么 g^i \bmod P 的结果两两不同,且有1<g<P, 0<i<P,归根到底就是 g^x \equiv 1 \pmod P 当且仅当

用原根代替\mathrm{FFT}中的单位根\omega_{n}^{k}\mathrm{NTT}常见模数998244353的原根是3

基础代码: ```cpp #include <iostream> #include <cstring> #include <algorithm> using namespace std; typedef long long ll; const int N=3000010; const int mod=998244353,G=3,niG=332748118; int n,m; int a[N],b[N]; int bit,tot,rev[N]; int qmi(ll a,int b,int p) { ll res=1; while (b) { if (b&1) res=res*a%p; a=a*a%p; b>>=1; } return res; } void NTT(int a[],int dir) { for (int i=0;i<tot;i++) { if (i<rev[i]) { swap(a[i],a[rev[i]]); } } for (int mid=1;mid<tot;mid<<=1) { int wn1=qmi((dir==1 ? G : niG),(mod-1)/(mid<<1),mod); for (int i=0;i<tot;i+=(mid<<1)) { int wnk=1; for (int j=0;j<mid;j++,wnk=(ll)wnk*wn1%mod) { int x=a[i+j],y=(ll)wnk*a[i+j+mid]%mod; a[i+j]=(x+y)%mod; a[i+j+mid]=(x-y+mod)%mod; } } } } int main() { cin >> n >> m; for (int i=0;i<=n;i++) cin >> a[i]; for (int i=0;i<=m;i++) cin >> b[i]; while ((1<<bit)<n+m+1) bit++; tot=(1<<bit); for (int s=0;s<tot;s++) rev[s]=(rev[s>>1]>>1)|((s&1)<<(bit-1)); NTT(a,1),NTT(b,1); for (int i=0;i<tot;i++) a[i]=(ll)a[i]*b[i]%mod; NTT(a,-1); int nitot=qmi(tot,mod-2,mod); for (int i=0;i<=n+m;i++) { int res=(ll)a[i]*nitot%mod; cout << res << " "; } cout << "\n"; return 0; } ``` ### 多项式求逆 满足$F(x)*G(x)=1$的多项式$F(x)$和$G(x)$互为乘法逆元 考虑倍增求解在模$x^n$意义下的多项式的逆元 在模$x^1$的意义下的多项式逆元相当于$x^0$的项的乘法逆元,可以直接求解,作为边界 在模$x^n$的意义下的乘法逆元可以用模$x^{\frac{n}{2}}$的意义下的结果计算得到。 假设 $$ F(x)*G^\prime(x)\equiv 1 \pmod {x^{\frac{n}{2}}} \\ F(x)*G(x)\equiv 1 \pmod {x^n} $$ 可以发现此时有 $$ F(x)*G(x)\equiv 1 \pmod {x^{\frac{n}{2}}} \\ F(x)*(G(x)-G^\prime(x))\equiv 0 \pmod {x^{\frac{n}{2}}} \\ G(x)-G^\prime(x)\equiv 0 \pmod {x^{\frac{n}{2}}} \\ (G(x)-G^\prime(x))^2\equiv 0 \pmod {x^{n}} \\ G(x)^2-2G(x)G^\prime(x)+{G^\prime(x)}^2\equiv 0 \pmod {x^{n}} \\ F(x)(G(x)^2-2G(x)G^\prime(x)+{G^\prime(x)}^2)\equiv 0 \pmod {x^{n}} \\ F(x)G(x)^2-2F(x)G(x)G^\prime(x)+F(x){G^\prime(x)}^2\equiv 0 \pmod {x^{n}} \\ G(x)-2G^\prime(x)+F(x){G^\prime(x)}^2 \equiv 0\pmod {x^n} \\ G(x)\equiv 2G^\prime(x)-F(x){G^\prime(x)}^2 \pmod {x^n} $$ 可以$\mathrm{NTT}$计算,每次$n$跳两倍,$O(n\log n)

封装:

namespace Poly
{
    typedef vector<int> poly;
    const int mod=998244353,G=3,niG=332748118;
    int bit,tot,rev[N];

    void print(poly a)
    {
        for (int i=0;i<a.size();i++) cout << a[i] << " ";
        cout << "\n";
    }

    void NTT(poly &a,int dir)
    {
        for (int i=0;i<tot;i++)
        {
            if (i<rev[i])
            {
                swap(a[i],a[rev[i]]);
            }
        }
        for (int mid=1;mid<tot;mid<<=1)
        {
            int wn1=qmi((dir==1 ? G : niG),(mod-1)/(mid<<1),mod);
            for (int i=0;i<tot;i+=(mid<<1))
            {
                int wnk=1;
                for (int j=0;j<mid;j++,wnk=(ll)wnk*wn1%mod)
                {
                    int x=a[i+j],y=(ll)wnk*a[i+j+mid]%mod;
                    a[i+j]=(x+y)%mod;
                    a[i+j+mid]=((ll)x-y+mod)%mod;
                }
            }
        }
    }

    void DFT(poly &a) 
    {
        NTT(a,1);
    }

    void IDFT(poly &a)
    {
        NTT(a,-1);
        int nitot=qmi(tot,mod-2,mod);
        for (int i=0;i<a.size();i++) a[i]=(ll)a[i]*nitot%mod;
    }

    poly operator + (poly a,poly b)
    {
        int len=max(a.size(),b.size());
        a.resize(len),b.resize(len);
        for (int i=0;i<len;i++) a[i]=(a[i]+b[i])%mod;
        return a;
    }

    poly operator - (poly a,poly b)
    {
        int len=max(a.size(),b.size());
        a.resize(len),b.resize(len);
        for (int i=0;i<len;i++) a[i]=((ll)a[i]-b[i]+mod)%mod;
        return a;
    }

    poly operator * (poly a,int b)
    {
        for (int i=0;i<a.size();i++) a[i]=(ll)a[i]*b%mod;
        return a;
    }

    poly operator * (poly a,poly b)
    {
        int n=a.size()-1,m=b.size()-1;
        bit=0;
        while ((1<<bit)<n+m+1) bit++;
        tot=(1<<bit);
        for (int i=0;i<tot;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));

        a.resize(tot),b.resize(tot);

        DFT(a),DFT(b);
        for (int i=0;i<tot;i++) a[i]=(ll)a[i]*b[i]%mod;
        IDFT(a);

        a.resize(n+m+1);

        return a;
    }

    poly Inv(poly F,int n)
    {
        poly G2,G;
        G2.resize(1);
        G2[0]=qmi(F[0],mod-2,mod);
        for (int k=2;;k<<=1)
        {
            // mod x^k
            // G(x)=2G2(x)-F(x)G2(x)G2(x) (mod x^n)

            poly t1=G2*G2,t2=F;
            t1.resize(k),t2.resize(k);
            poly t=t1*t2;
            t.resize(k);
            G=G2*2-t;
            G.resize(k);
            G2=G;
            if (k>=n) break;
        }
        G.resize(n+1);

        return G;
    }
};

原题: P4841

正解:

$$ f_i=g_i-\sum_{j=1}^{i-1}\binom{i-1}{j-1}f_jg_{i-j} \\ g_i=f_i+\sum_{j=1}^{i-1}\binom{i-1}{j-1}f_jg_{i-j} \\ g_i=\sum_{j=1}^{i}\binom{i-1}{j-1}f_jg_{i-j} \\ g_i=\sum_{j=1}^{i}\frac{(i-1)!}{(j-1)!(i-j)!}f_jg_{i-j} \\ \frac{g_i}{(i-1)!}=\sum_{j=1}^{i}\frac{f_jg_{i-j}}{(j-1)!(i-j)!} \\ \frac{g_i}{(i-1)!}=\sum_{j=1}^{i}\frac{f_j}{(j-1)!}\frac{g_{i-j}}{(i-j)!} \\ A(x)=\sum_{i=1}\frac{f_i}{(i-1)!}x^i \\ B(x)=\sum_{i=1}\frac{g_i}{(i-1)!}x^i \\ C(x)=\sum_{i=0}\frac{g_i}{i!}x^i \\ B=A*C \\ A=B*C^{-1} $$ 有$g_i=2^{\frac{n(n-1)}{2}}$,直接预处理$g_i$和阶乘逆元,用多项式乘法和多项式求逆计算,复杂度$O(n\log n)