多项式

多项式杂谈(更新中)

前言

最近在学习多项式,稍微开个坑吧。

一些比较好的学习资料、模板:

我个人习惯写指针,在我的模板中基本没有 vector,可能是我觉得 vector 常数大吧。

本文以存板子为主,一些边界条件可以参照洛谷模板题的输入格式。

快速傅里叶变换 FFT

题目大意:给定两个多项式 \(f(x),g(x)\),求 \(f(x) * g(x)\)。即求:\(h[i]=\sum f[j]\times g[i-j]\)。

模板题 学习资料

这题显然有一个 \(O(n^2)\) 的暴力,不过我们有一个 \(O(n\log n)\) 的优(du)(liu)算法,常数大就是了。

一些前置知识:复数,单位根的性质。

我们知道 \(n\) 个点可以确定一个 \(n-1\) 次多项式,即用 \(n\) 个点 \((x_i,y_i)\) 来表示一个多项式。我们称之为点值表示法。而两个点值多项式乘起来只需要将对应点值相乘即可,即 \((x_i,y1_i), (x_i,y2_i)\to (x_i,y1_i\times y2_i)\),时间复杂度是 \(O(n)\) 的。

注意 \(n\) 次多项式和 \(m\) 次多项式相乘要选 \(n+m+1\) 个点值(多项式的系数是 \(n+m\) 次的)。

接下来我们的问题就是如何将系数多项式和点值多项式互相转化。我们将系数 \(\to\) 点值的过程称为 DFT(离散傅里叶变换),将点值 \(\to\) 系数的过程称为 IDFT (离散傅里叶逆变换)。

如果我们暴力带入点值,时间复杂度还是 \(O(n^2)\) 的,不可接受。但是我们可以通过带入一些特殊的点值来优化这个过程。下面不妨假设 \(n\) 为 \(2\) 的次幂。设一个多项式 \(A(x)\)。

\[ {\displaystyle} A(x)=\sum _{i=0} ^ {n-1} a_ix^i \]

将 \(x^i\) 按 \(i\) 的奇偶性分成两部分,并将后一部分的 \(x\) 提出来,有:

\[ {\displaystyle} A(x)=\sum _{i=0} ^{\lfloor \frac n2 \rfloor-1} a_{2i}x^{2i}+x\sum _{i=0} ^{\lfloor \frac n2 \rfloor-1} a_{2i+1}x^{2i} \]

再设两个多项式 \(A_1(x),A_2(x)\)。

\[A_1(x)=a_0x^0+a_2x^2+\cdots \\ A_2(x)=a_1x^0+a_3x^2+\cdots \\ A(x)=A_1(x^2)+xA_2(x^2) \]

设 \(k< \frac n2\),将 \(\omega_n^k\) 和 \(\omega_n^{k+n/2}\) 带入。

\[A(\omega_n^k)=A_1((\omega_n^k)^2)+\omega_n^kA_2((\omega_n^k)^2) \\ =A_1(\omega_{n/2}^k)+\omega_n^kA_2(\omega_{n/2}^k) \\ A(\omega_n^{k+n/2})=A_1((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}A_2((\omega_n^{k+n/2})^2) \\ =A_1(\omega_{n/2}^k)-\omega_n^kA_2(\omega_{n/2}^k) \]

注意到上下两个式子只有符号不同,也就是说,我们带入 \(\omega_n^k\) 求出 \(A_1(\omega_n^k)\) 和 \(A_2(\omega _n^k)\) 就能得到两个点值。

这样我们就可以递归求解了,时间复杂度 \(T(n)=2T(\frac n2)+O(n)=O(n\log n)\)。

接下来解决 IDFT

我们得到了点值 \(s[k]=\sum \limits _ {i=0} ^ {n-1} f[i]\times (\omega _ n ^ k) ^ i\)。要求的多项式 \(n\times p[k]=\sum \limits _ {i=0} ^{n-1} s[i]\times (\omega_n^{n-k})^i\)。证明带入即可。这样 IDFT 就和 DFT 一样了,唯一的区别就是带入了逆单位根。

这样我们就可以写出递归版的 FFT 了,但是常数极大,可能模板题都过不了。我们需要些非递归版的 FFT

但是在 FFT 的时候需要将位置分奇偶,这样和原来的位置就不一样了。我们可以先将原序列变换好,放在最终的位置上,再一步一步向上合并。

注意到每个数最后的位置为其二进制翻转后的得到的位置。

struct Complex {
    double x,y;
    inline Complex(double a=0,double b=0) {x=a,y=b;}
    inline Complex operator + (const Complex a) const {return Complex(x+a.x,y+a.y);}
    inline Complex operator - (const Complex a) const {return Complex(x-a.x,y-a.y);}
    inline Complex operator * (const Complex a) const {return Complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}b[N],a[N];
int rev[N],f[N],g[N];
inline void FFT(Complex *a,int len,int f) {
    int k=0;
    while ((1<<k)<len) ++k;
    for (int i=0;i<len;i++) {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
        if (i<rev[i]) swap(a[i],a[rev[i]]);
    }
    for (int l=2;l<=len;l<<=1) {
        Complex wn(cos(2*pi*f/l),sin(2*pi*f/l));
        for (int j=0;j<len;j+=l) {
            Complex w(1.0,0.0);
            for (int k=j;k<j+l/2;k++) {
                Complex t=w*a[k+l/2];
                a[k+l/2]=a[k]-t,a[k]=a[k]+t,w=w*wn;
            }
        }
    }
}
inline void DFT(Complex *a,int len) {
    FFT(a,len,1);
}
inline void IDFT(Complex *a,int len) {
    FFT(a,len,-1);
    for (int i=0;i<=len;i++) a[i].x/=(double)len;
}
inline void Mul(Complex *a,Complex *b,int n,int m) {
    int len=1; while (len<=(n+m)) len<<=1;
    for (int i=n;i<len;i++) a[i]=Complex(0,0);
    for (int i=m;i<len;i++) b[i]=Complex(0,0);
    DFT(a,len),DFT(b,len);
    for (int i=0;i<=len;i++) a[i]=a[i]*b[i];
    IDFT(a,len);
}

快速数论变换 NTT

FFT 要算一大堆三角函数常数巨大,还可能会出现精度问题导致 WA。这个时候需要 NTT

题目大意:给定两个多项式 \(f(x),g(x)\),求 \(f(x) * g(x)\)。即求:\(h[i]=\sum f[j]\times g[i-j]\)。对 \(998244353\) 取模。

这道模板题中最终系数保证不会大于 \(998244353\),所以可以视为对 \(998244353\) 取模。

模板题 学习资料

前置知识:原根。

FFT 之所以快是因为 \(\omega\) 有优秀的性质,实际上原根也有同样的性质。相当于将 \(\omega_n^1\) 变成 \(g^{\frac {p-1}n}\)。其它部分都差不多。

但是需要注意:设模数为 \(p=r\times 2^k+1\)。NTT 只能处理 \(n\leq 2^k\) 的数据,比如模数为 \(998244353=119\times 2^{23}+1\),此时原根为 \(3\)。但当模数比较毒瘤,比如 \(10^9+7\),此时 \(k=1\) ,NTT 不能再使用。具体可以参照这篇博客

int rev[N],a[N],b[N],c[N];
inline void NTT(int *a,int len,int f) {
    int k=0;
    while ((1<<k)<len) ++k;
    for (int i=0;i<len;i++) {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
        if (i<rev[i]) swap(a[i],a[rev[i]]);
    }
    for (int l=2;l<=len;l<<=1) {
        int m=l>>1,w=Pow(G,(mod-1)/l);
        if (f) w=Pow(w,mod-2);
        for (int i=0;i<len;i+=l) {
            int wk=1;
            for (int j=0;j<m;j++,wk=1LL*wk*w%mod) {
                int p=a[i+j],q=1LL*wk*a[i+j+m]%mod;
                a[i+j]=(p+q)%mod,a[i+j+m]=(p-q+mod)%mod;
            }
        }
    }
}
inline void DFT(int *a,int len) {
    NTT(a,len,false);
}
inline void IDFT(int *a,int len) {
    NTT(a,len,true); int invn=Pow(len,mod-2);
    for (int i=0;i<len;i++) a[i]=1LL*a[i]*invn%mod;
}

多项式求逆

题目大意,给定多项式 \(f(x)\) ,求一个多项式 \(g(x)\) ,满足 \(f(x)*g(x)\equiv 1\pmod {x^n}\)。

模板题

主要的思想是倍增。假设我们已知 \(h(x)*f(x)\equiv 1 \pmod {x^{n/2}}\)。

有 \(g(x)-h(x)\equiv 1 \pmod {x^{n/2}}\)。

两边平方,有: \(g^2(x)-2g(x)h(x)+h^2(x)\equiv 1 \pmod {x^{n}}\)。

两边同乘 \(f(x)\),有:\(g(x)-2h(x)+h^2(x)f(x)\equiv 1 \pmod {x^n}\)。

所以有:\(g(x)=2h(x)-h^2(x)f(x)+1\pmod {x^n}\)。

倍增求解 \(g(x)\) 即可。

inline void PolyInv(int n,int *a,int *b) {
    if (n==1) return b[0]=Pow(a[0],mod-2),void();
    PolyInv((n+1)>>1,a,b); int len=1;
    while (len<(n<<1)) len<<=1;
    for (int i=0;i<n;i++) c[i]=a[i];
    for (int i=n;i<len;i++) c[i]=0;
    DFT(c,len),DFT(b,len);
    for (int i=0;i<len;i++)
        b[i]=1LL*(2-1LL*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
    IDFT(b,len);
    for (int i=n;i<len;i++) b[i]=0;
}

多项式除法&取模

多项式对数函数

多项式指数函数

多项式快速幂

多项式多点求值

多项式快速插值

一些卡常技巧和提示

多项式题一般比较难写,而且常数大。下面给出一些卡常技巧和提示。

  • Mul(a,b) 的时候最好将 ab 先放到临时变量中。

  • 注意变量不要重名,static 可以很好的解决这个问题,并优化内存,但是可能会增大常数。

  • 在做一些多项式操作前一定要想清楚边界范围,并记得清空数组。这点可能使用 vector 会比较方便。

  • 在做一些多项式操作前一定要记得清空数组。

  • 使用取模优化

  • 加一些 inlineregister

  • static 好像会变慢?

  • 使用一些操作减少 FFT 次数,比如三次变两次的技巧。

  • 预处理可能要用到的所有原根次幂,这样在求原根次幂的时候回访问一段连续内存加快速度。具体可以如下:

    int GPow[2][20][N];
    inline void InitG() {
        for (int p=1;p<=19;p++) {
            int buf1=Pow(G,(Mod-1)/(1<<p));
            int buf0=Pow(invG,(Mod-1)/(1<<p));
            GPow[1][p][0]=GPow[0][p][0]=1;
            for (int i=1;i<(1<<p);i++)
                GPow[1][p][i]=1LL*GPow[1][p][i-1]*buf1%Mod,
                GPow[0][p][i]=1LL*GPow[0][p][i-1]*buf0%Mod;
        }
    }
    

    同时将 NTT 改成这样

    inline void NTT(int *a,int len,int f) {
        int k=0;
        while ((1<<k)<len) ++k;
        for (Re int i=0;i<len;++i)
            if (i<rev[i]) swap(a[i],a[rev[i]]);
        for (Re int l=2,cnt=1;l<=len;l<<=1,cnt++) {
            int m=l>>1;
            for (Re int i=0;i<len;i+=l) {
                int *buf=GPow[f^1][cnt];
                for (Re int j=0;j<m;j++,buf++) {
                    int p=a[i+j],q=1LL*(*buf)*a[i+j+m]%Mod;
                    a[i+j]=chk(p+q),a[i+j+m]=chk(p-q+Mod);
                }
            }
        }
    }
    

结语

咕咕咕

上一篇:(七)根轨迹①


下一篇:(七)二阶系统对初始条件的动态响应