多项式

Posted 2020-12-06 lenstar

多项式杂谈（更新中）

前言

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

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

• 多项式大杂烩 - Froggy 珂学家的博客
• 【多项式n连】各种多项式模板（从入门到入土）

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

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

快速傅里叶变换 FFT

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

模板题 学习资料

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

一些前置知识：复数，单位根的性质。

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

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

接下来我们的问题就是如何将系数多项式和点值多项式互相转化。我们将系数 ( o) 点值的过程称为 DFT（离散傅里叶变换），将点值 ( o) 系数的过程称为 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 floor-1} a_{2i}x^{2i}+xsum _{i=0} ^{lfloor frac n2 floor-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(nlog n))。

接下来解决 IDFT 。

我们得到了点值 (s[k]=sum limits _ {i=0} ^ {n-1} f[i] imes (omega _ n ^ k) ^ i)。要求的多项式 (n imes p[k]=sum limits _ {i=0} ^{n-1} s[i] imes (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] imes g[i-j])。对 (998244353) 取模。

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

模板题 学习资料

前置知识：原根。

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

但是需要注意：设模数为 (p=r imes 2^k+1)。NTT 只能处理 (nleq 2^k) 的数据，比如模数为 (998244353=119 imes 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 1pmod {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)+1pmod {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) 的时候最好将 a 和 b 先放到临时变量中。

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

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

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

• 使用取模优化

• 加一些 inline 和 register。

• 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);
              }
          }
      }
  }
  

结语

咕咕咕