-1.原根的引入

FFT 对单位根的依赖，决定该算法必须采用浮点数，进而引发精度问题。

糟糕的是，数学家们已经证明，在复数域 $\mathbb C$ 中，单位根是唯一一类满足要求的数。

大部分的计数题都是在模意义下完成的，我们自然希望为 FFT 找一个模意义下的替代品。

考虑引入模意义同余系 $F_p$ 中单位根的类似物 : 原根。

可能有帮助的文章 : 原根&离散对数相关

我们先罗列一下 FFT 中用到的 $ω_n$ 的所有性质。

① $ω_n^k=(ω_n^1)^k$

这正是单位根定义的一部分。

由此，我们只需要找出一个满足要求的 $ω_n^1$ 即可。
② $ω_n^{0\sim(n-1)}$ 必须互不相同。

否则点值就会有重复，插值就会不正确。
③ $ω_n^k=ω_n^{k\bmod n}$

可以导出对称引理 : $\omega_{n}^{k+n/2}=-\omega_{n}^k$

当 $n$ 是偶数的时候，通过对称引理即能够得出求和引理。

综合①②③，等价于 : $w_n^1$ 在模意义下的阶恰为 $n$ 。
④ $\omega_{2n}^{2k}=\omega_{n}^{k}$ : 折半引理。

等价于 $(\omega_{2n})^2=\omega_n$ 。

我们来看看原根的定义。

众所周知，对于素数 $p$ ，其模意义同余系 $F_p$ 中恰有 $1,2,3,...p-1$ 这 $p-1$ 个数。

对于 $g$ ，若 $g$ 的阶达到 $p-1$ 的上界，则称 $g$ 为原根。

显然， $p-1$ 并不一定等于 $n$ ，所以，原根本身并不能直接用作单位根。

但是，能够发现， $g^k$ 的阶数恰为 $(p-1)/\gcd(k,p-1)$ ，这个数仍然是 $p-1$ 的约数。

所以，当 $n\!\not|\ (p-1)$ 时，找不到阶恰为 $n$ 的数。

当 $n|(p-1)$ 时， $g^{(p-1)/n}$ 的阶则恰为 $(p-1)/\gcd((p-1)/n,p-1)=n$ 。

这说明， $g^{(p-1)/n}$ 满足了要求①②③。

④ : $(\omega_{2n})^2=(g^{(p-1)/2n})^{2}=g^{(p-1)/2n*2}=g^{(p-1)/n}=(\omega_{n}^1)^{k}=\omega_n$ ，得证。

所以， $g^{(p-1)/n}$ 符合我们的所有需求。

前文提到过，当且仅当 $n|(p-1)$ 才能构造出单位根。事实上，分治算法中的 $n$ 往往是 $2$ 的幂，我们只需要让 $p-1$ 包含大量因子 $2$ 即可。

比如 $998244353=2^{23}*7*17+1$ 即为一个很有代表性的模数。

熟练的选手会记下一些常用模数的原根。

懒人福利 : 安利大佬的原根表。

至于如何寻找原根，可见上面拉链的那篇文章。

但是，一定非要原根不可吗?

摆脱原根的方法来自 Uoj : hly1204

假如我们已经知道 $p-1=s*2^r$ 且 $s$ 为奇数。

随机出一个二次非剩余 $v$ ，若 $g$ 为原根且有 $v=g^k$ ，那么 $k$ 一定是个奇数。

那么有 ${\rm ord}(v^s)={\rm ord}(g^{ks})=\dfrac{s*2^r}{\gcd(s*2^r,sk)}=2^r$。

则有 $\omega_{n}=(v^s)^{2^r/n}$ 。

0.NTT

好的，我们找到单位根的替代品了，现在我们来改代码。

上次讲的 FFT 的最终版本(没有“三次变两次优化”):

Old评测记录

#include<algorithm>
#include<cstdio>
#include<cmath>
#define Maxn 1350000
using namespace std;
const double Pi=acos(-1);
int n,m;
struct CP
{
  CP (double xx=0,double yy=0){x=xx,y=yy;}
  double x,y;
  CP operator + (CP const &B) const
  {return CP(x+B.x,y+B.y);}
  CP operator - (CP const &B) const
  {return CP(x-B.x,y-B.y);}
  CP operator * (CP const &B) const
  {return CP(x*B.x-y*B.y,x*B.y+y*B.x);}
}f[Maxn<<1],p[Maxn<<1];
int tr[Maxn<<1];
void fft(CP *f,bool flag)
{
  for (int i=0;i<n;i++)
    if (i<tr[i])swap(f[i],f[tr[i]]);
  for(int p=2;p<=n;p<<=1){
    int len=p>>1;
    CP tG(cos(2*Pi/p),sin(2*Pi/p));
    if(!flag)tG.y*=-1;
    for(int k=0;k<n;k+=p){
      CP buf(1,0);
      for(int l=k;l<k+len;l++){
        CP tt=buf*f[len+l];
        f[len+l]=f[l]-tt;
        f[l]=f[l]+tt;
        buf=buf*tG;
      }
    }
  }
}
int main()
{
  scanf("%d%d",&n,&m);
  for (int i=0;i<=n;i++)scanf("%lf",&f[i].x);
  for (int i=0;i<=m;i++)scanf("%lf",&p[i].x);
  for(m+=n,n=1;n<=m;n<<=1);
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
  fft(f,1);fft(p,1);//DFT
  for(int i=0;i<n;++i)f[i]=f[i]*p[i];
  fft(f,0);
  for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.49));
  return 0;
}

摇身一变!NTT出现!

还是原来的配方,还是熟悉的味道。

Code

#include<algorithm>
#include<cstdio>
#define ll long long
#define mod 998244353
#define G 3
#define Maxn 1350000
using namespace std;
inline int read()
{
  register char ch=0;
  while(ch<48||ch>57)ch=getchar();
  return ch-'0';
}
ll powM(ll a,ll t=mod-2)
{
  ll ans=1;
  while(t){
  	if(t&1)ans=ans*a%mod;
  	a=a*a%mod;t>>=1;
  }return ans;
}
int n,m,tr[Maxn<<1];
ll f[Maxn<<1],g[Maxn<<1],invn;
const ll invG=powM(G);
void NTT(ll *f,bool op)
{
  for (int i=0;i<n;i++)
    if (i<tr[i])swap(f[i],f[tr[i]]);
  for(int p=2;p<=n;p<<=1){
    int len=p>>1,
        tG=powM(op ? G:invG,(mod-1)/p);
    for(int k=0;k<n;k+=p){
      ll buf=1;
      for(int l=k;l<k+len;l++){
        int tt=buf*f[len+l]%mod;
        f[len+l]=f[l]-tt;
        if (f[len+l]<0)f[len+l]+=mod;
        f[l]=f[l]+tt;
        if (f[l]>mod)f[l]-=mod;
        buf=buf*tG%mod;
      }
    }
  }
}
int main()
{
  scanf("%d%d",&n,&m);n++;m++;
  for (int i=0;i<n;i++)f[i]=read();
  for (int i=0;i<m;i++)g[i]=read();
  for(m+=n,n=1;n<m;n<<=1);
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
  NTT(f,1);NTT(g,1);
  for(int i=0;i<n;++i)f[i]=f[i]*g[i]%mod;
  NTT(f,0);invn=powM(n);
  for(int i=0;i<m-1;++i)
    printf("%d ",(int)(f[i]*invn%mod));
  return 0;
}

评测记录

未经毒瘤卡常时,单个 NTT 变换的时间为 FFT 的 75%~80%.

其中这一句话改变得很大:

ll tG=powM(op==1 ? G:invG,(mod-1)/p);

即 $\omega_p^1=g^{(mod-1)/p}$ ,再也不用三角函数啦!

一些额外的技巧:

预处理单位根，这样可以省去部分乘法操作，并访问一段连续内存。
注意到只会有加减法操作，可以使用 ull 存储，能耐受 18*mod*mod 的范围，在常规范围下可以不取模。

模板题范围较大,在中间取一次模即可。

在不同的题目中，可以加速 15%~40% 不等。

#include<algorithm>
#include<cstdio>
#define ll long long
#define ull unsigned long long
using namespace std;
const int _G=3,mod=998244353,Maxn=1050000;
inline int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
ll powM(ll a,int t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1];
void NTT(int *g,bool op,int n)
{
  static ull f[Maxn<<1],w[Maxn<<1]={1};
  for (int i=0;i<n;i++)f[i]=g[tr[i]];
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }
    if (l==(1<<17))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
int n,m,F[Maxn<<1],G[Maxn<<1];
int main()
{
  n=read()+1;m=read()+1;
  for (int i=0;i<n;i++)F[i]=read();
  for (int i=0;i<m;i++)G[i]=read();
  int len=1;for (;len<n+m;len<<=1);
  for(int i=0;i<len;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?len>>1:0);
  NTT(F,1,len);NTT(G,1,len);
  for(int i=0;i<len;++i)F[i]=1ll*F[i]*G[i]%mod;
  NTT(F,0,len);
  for (int i=0;i<n+m-1;i++)
    printf("%d ",(int)F[i]);
  return 0;
}

$\frac{1}{2}$ .注意事项

关于封装

多项式模板丰富之后，如何封装成为令人头痛的大难题。

前面讲过，利用 FFT(NTT) 变换的线性性可以减少运算。

这样就需要适度封装以利用性质卡常数，下面给出比较智能的封装:

IO (仅模板)

inline int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
inline void print(ll *f,int len){
  for (int i=0;i<len;i++)
    printf("%lld ",f[i]);
  puts("");
}

多项式基本操作

#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
const int _G=3,mod=998244353,Maxn=;

ll powM(ll a,ll t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn<<1],w[Maxn<<1]={1};
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }
    if (l==(1<<10))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
void times(int *f,int *g,int len,int lim)
{
  static int sav[Maxn<<1];
  int n=1;for(n;n<len+len;n<<=1);
  clr(sav,n);cpy(sav,g,n);
  NTT(f,1,n);NTT(sav,1,n);
  px(f,sav,n);NTT(f,0,n);
  clr(f+lim,n-lim);clr(sav,n);
}

熟练的选手能够在 8~15min 内默完，保险起见建议拿 times 跟暴力卷积拍一下。

这样我们就能够肆无忌惮的使用 NTT 和一系列缩写了。

小心地分析清零带来的影响，宁愿多清不能漏网，而且这种问题一般很难检查和拍出来。

善用宏可以有效避免重名问题。

关于运算合法性/理解方式

本文所讨论的多项式均为形式幂级数，但限于篇幅，并没有给出其严格定义和应用。

它只是辅助分析的数学工具，所谓多项式中的“未知数”并没有实际意义，只是一个代数对象。

如果想了解更多相关内容，可见 rqy : 浅谈 OI 中常用的一些生成函数运算的合法与正确性。

简化地讲，我们将以多项式加法和加法卷积来定义所有运算。

在大多数题目中，我们只对多项式的前若干项感兴趣，此时我们为所有运算设定一个次数上界，即 $\pmod{x^n}$ 。

由于多项式加法和乘法转移的单向性（只会从低次向高次贡献），不难发现下列两条性质：

$$\big(F(x)\bmod{x^n}+(G(x)\bmod{x^n}\big)\bmod{x^n}=\big(F(x)+G(x)\big)\bmod{x^n}$$$$\big(F(x)\bmod{x^n}\big)*\big(G(x)\bmod{x^n}\big)\bmod{x^n}=\big(F(x)*G(x)\big)\bmod{x^n}$$

前面提到过，所有运算都是利用加法和乘法定义的，所以这能说明，我们可以在每一步运算中都保持次数界，利用其去除无用的高次项。

这样就能避免对无穷幂级数的处理，以保障复杂度。

1.多项式四则运算

终于从乘法向多项式全家桶迈进了!

约定： $[x^n]F(x)$ 和 $F[n]$ 均表示 $F(x)$ 的 $n$ 次项系数。

也就是说， $F(x)$ 可写作 $\sum\limits_{i=0}F[i]x^i$ 。

多项式的加减定义为：

F(x)+G(x)=\sum\limits_{i=0}\big(F[i]+G[i]\big)x^i

F(x)-G(x)=\sum\limits_{i=0}\big(F[i]-G[i]\big)x^i

其实就是简单地将对应项相加减， $O(n)$ 即可完成计算。

多项式乘法（加法卷积）定义为：

$$F(x)*G(x)=\sum\limits_{i=0}\sum\limits_{j=0}F[i]G[i]x^{i+j}=\sum\limits_{i=0}x^i\sum\limits_{k=0}^iF[k]G[i-k]$$

使用 $NTT$ 在 $O(n\log n)$ 内完成计算。

目前为止，加法，减法，乘法都和我们平时多项式运算的经验保持一致。

接下来就是多项式乘法逆元，定义为：

满足 $F(x)*G(x)=1$ 时，称 $F(x),G(x)$ 互为乘法逆元。

类比同余系下数字的逆元，多项式逆元可以把除法转换成乘法。

模板：P4238 【模板】多项式乘法逆

怎么求一个多项式的乘法逆元呢？

方法一：递推

由 $F(x)G(x)=1$ ，可以得到:

$\sum\limits_{i=0}^nF[i]G[n-i]=[n=0]$ , 下面假定 $n>0$

$\sum\limits_{i=0}^{n-1}F[i]G[n-i]+F[n]G[0]=0$

$F[n]=-\frac{1}{G[0]}\sum\limits_{i=0}^{n-1}F[i]G[n-i]$

可以 $O(n)$ 递推出末项，计算整个乘法逆元的复杂度即为 $O(n^2)$ 。

这个递推算法并没有用到 NTT，所以不可能达到 $O(n\log n)$ 的理想复杂度。
方法二：倍增

我们考虑采用倍增的思想，不断扩大次数上界。

当界为 $\pmod {x^1}$ ，即多项式只有常数项的时候，直接求数字逆元即可满足条件。

假设我们已经得到了 $R_*(x)=F(x)^{-1}(mod\ x^{n/2})$ ，即逆元的前 $n/2$ 位。

设 $R(x)=F(x)^{-1}(mod\ x^n)$ ，即逆元的前 $n$ 位。

明显有 $R(x)=R_*(x)(mod\ x^{n/2})$ ，做差得到 $R(x)-R_*(x)=0(mod\ x^{n/2})$

如何把 $(mod\ x^{n/2})$ 提升到 $(mod\ x^n)$ 呢? 平方。

得 $(R(x)-R_*(x))^2=0(mod\ x^n)$

展开得到 $R(x)^2-2R(x)R_*(x)+R_*(x)^2=0(mod\ x^n)$

等式两边同乘 $F(x)$ 得到 $R(x)-2R_*(x)+R_*(x)^2F(x)=0(mod\ x^n)$

移项得到 $R(x)=2R_*(x)-R_*(x)^2F(x)(mod\ x^n)$

根据这个即可从 $R_*(x)$ 求得 $R(x)$ ,实现次数倍增。

复杂度为 $T(n)=T(n/2)+O(n\log n)=O(n\log n)$ ，下面所有倍增法的复杂度类似。

Code : (配合前面的封装食用)

void invp(int *f,int m)
{
  int n;for (n=1;n<m;n<<=1);
  static int w[Maxn<<1],r[Maxn<<1],sav[Maxn<<1];
  w[0]=powM(f[0]);
  for (int len=2;len<=n;len<<=1){
  	for (int i=0;i<(len>>1);i++)
      r[i]=(w[i]<<1)%mod;
    cpy(sav,f,len);
    NTT(w,1,len<<1);px(w,w,len<<1);
    NTT(sav,1,len<<1);px(w,sav,len<<1);
  	NTT(w,0,len<<1);clr(w+len,len);
  	for (int i=0;i<len;i++)
      w[i]=(r[i]-w[i]+mod)%mod;
  }cpy(f,w,m);clr(sav,n+n);clr(w,n+n);clr(r,n+n);
}
int n,F[Maxn<<1];
int main()
{
  n=read();
  for (int i=0;i<n;i++)F[i]=read();
  invp(F,n);print(F,n);
  return 0;
}

注意,由于中间 NTT 是两倍长度,最后清空数组要清空两倍的 n。

检查的方法 : 随机一个多项式求逆两次看看是不是自己。

评测记录

卡常后的倍增

众所周知 $\rm NTT$ 其实是在计算长度为 $2^k$ 的循环卷积，我们可以利用这个性质来简化计算。

观察倍增公式 $R(x)=2R_*(x)-R_*(x)^2F(x)(mod\ x^n)$ ，我们需要使用乘法的项仅有 $R_*(x)^2F(x)$ 。

我们最终只需要得到答案的后 $n/2$ 项。

先计算 $F(x)*R_*(x)$ ，两者的次数分别为 $n,n/2$ 。

这部分结果的前 $n/2$ 项都为 $0$ ，所以，循环卷积溢出只会破坏前面的 $0$ 。

再乘上一个 $R_*(x)$ ，此时循环卷积溢出指挥破坏前 $n/2$ 项，正好是我们不需要的。

void invp(int *f,int m)
{
  int n;for (n=1;n<m;n<<=1);
  static int w[Maxn<<1],r[Maxn<<1],sav[Maxn<<1];
  w[0]=powM(f[0]);
  for (int len=2;len<=n;len<<=1){
  	for (int i=0;i<(len>>1);i++)r[i]=w[i];
    cpy(sav,f,len);NTT(sav,1,len);
    NTT(r,1,len);px(r,sav,len);
    NTT(r,0,len);clr(r,len>>1);
    cpy(sav,w,len);NTT(sav,1,len);
  	NTT(r,1,len);px(r,sav,len);
  	NTT(r,0,len);
  	for (int i=len>>1;i<len;i++)
      w[i]=(w[i]*2ll-r[i]+mod)%mod;
  }cpy(f,w,m);clr(sav,n);clr(w,n);clr(r,n);
}

评测记录 (可以加速两倍多)

2.多项式牛顿迭代(及求导/积分/复合)

继续学习之前，建议读者先自行了解简单微积分。

首先定义多项式的求导 :

F'(x)=\sum\limits_{i=0}F[i]ix^{i-1}

定义多项式的积分 :

F'(x)=C+\sum\limits_{i=0}F[i]x^{i+1}/i

这和经典的求导、积分是一致的。

Code:

void dao(int *f,int m){
  for (int i=1;i<m;i++)
    f[i-1]=1ll*f[i]*i%mod;
  f[m-1]=0;
}
int inv[Maxn];
void Init(int lim){
  inv[1]=1;
  for (int i=2;i<=lim;i++)
    inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
void jifen(int *f,int m){
  for (int i=m;i;i--)
    f[i]=1ll*f[i-1]*inv[i]%mod;
  f[0]=0;
}

定义多项式复合为：

F(G(x))=\sum\limits_{i=0}F[i]G(x)^i

将所有 $x$ 代换为 $G(x)$ 即可得到上式，利用整体法不难理解。

多项式牛顿迭代

用于解决下列问题：

已知函数 $G$ 且 $G(F(x))=0$ ，求多项式 $F\pmod {x^n}$ 。

实际应用中函数 $G$ 一般较简单，而且是个常量(可以手动数学分析)。

证明可在多项式计数杂谈中查看。

结论：$\Large{F(x)=F_*(x)-\frac{G(F^*(x))}{G'(F^*(x))}}\pmod{x^n}$

也就是说我们采用这个式子倍增地求解上述问题。

这个式子很牛逼，可以帮我们省下思维难度(肝)，一定好好记住。

需要注意的是， $G(F^*(x))$ 的前 $n/2$ 项均为 $0$ ，所以分母的计算精度就只需要达到 $\pmod {x^{n/2}}$ 级别。

比如说上面的 多项式求逆

根据 $B(x)*A(x)-1=0(mod\ x^n)$ ,这里 $A(x)$ 已知。

设 $G(x)$ 且 $G(B(x))=A(x)B(x)-1$ ，这里 $A(x)$ 是系数。

所以 $G'(B(x))=A(x)(mod\ x^n)$ (注意这个求导的主元是 $B(x)$ 而非 $x$ )

(下文默认 $B_*(x)$ 表示在 $(mod\ x^\frac{n}{2})$ 处的解)。

得 $B(x)=B_*(x)-\frac{G(B_*(x))}{G'(B_*(x))}=B_*(x)-\frac{1}{A(x)}(A(x)B_*(x)-1)$

注意 $\frac{1}{A(x)}$ 的精度只需要达到 $\pmod {x^{n/2}}$ ，所以直接使用 $B_*(x)$ 即可。

得 $B(x)=B_*(x)-B_*(x)(A(x)B_*(x)-1)=B_*(x)(2-A(x)B_*(x))$

这与我们上面采用平方法推出的式子是等价的。
多项式开方

P5205 【模板】多项式开根

得 $B(x)^2-A(x)=0$ ，此处 $A(x)$ 已知。

令 $G(B(x))=B(x)^2-A(x)$

则 $G'(B(x))=2B(x)$

根据牛顿迭代得到 : $B(x)=B^*(x)-\frac{G(B^*(x))}{G'(B^*(x))}$

$=B^*(x)-\frac{B^*(x)^2-A(x)}{2B^*(x)}=B^*(x)+\frac{A(x)-B^*(x)^2}{2B^*(x)}=\frac{A(x)+B^*(x)^2}{2B^*(x)}$

根据这个倍增即可。

多项式只有常数时，开方的结果是二次剩余，题目中保证常数为 $1$ ，那么直接令 $B_0=1$ 。

Code:

void sqrtp(int *f,int m)
{
  int n;for (n=1;n<m;n<<=1);
  static int b1[Maxn<<1],b2[Maxn<<1];
  b1[0]=1;
  for (int len=2;len<=n;len<<=1){
  	for (int i=0;i<(len>>1);i++)b2[i]=(b1[i]<<1)%mod;
  	invp(b2,len);
  	NTT(b1,1,len);px(b1,b1,len);NTT(b1,0,len);
  	for (int i=0;i<len;i++)b1[i]=(f[i]+b1[i])%mod;
    times(b1,b2,len,len);
  }cpy(f,b1,m);clr(b1,n+n);clr(b2,n+n);
}

评测记录。

3.多项式带余除法

P4512 【模板】多项式除法

如果保证没有余多项式（即整除），直接使用普通多项式除法就能得到商多项式。

我们考虑用一种奥妙重重的方式把余多项式去掉。

前面我们经常利用 $\pmod {x^?}$ 来消去某些项。很明显，这种操作只能 消去高次的项,而留下低次的项。

然而此处余多项式的次数低，考虑反转系数，将原来次数高的变成次数低的。

我们定义 $n$ 次多项式翻转操作 $F^T(x)=x^nF(x^{-1})$ 。

手玩一下就会发现 $F^T(x)$ 的系数是 $F(x)$ 的反转。

比如说 $F(x)=3x^2+2x+1$ , $\ F^T(x)=x^3F(x^{-1})=x^3(3x^{-2}+2x^{-1}+1)=x^3+2x^2+3$

首先有 $F(x)=Q(x)*G(x)+R(x)$ ,其中 $F(x)$ 和 $G(x)$ 已知。

换元得到 $F(x^{-1})=Q(x^{-1})*G(x^{-1})+R(x^{-1})$

两边同乘 $x^n$ : $x^nF(x^{-1})=x^nQ(x^{-1})*G(x^{-1})+x^nR(x^{-1})$

写成反转的形式，得 : $F^T(x)=Q^T(x)*G^T(x)+x^{n-m+1}R^T(x)$

此时我们机智地 $\bmod\ x^{n-m+1}$ ，进而得到 $F^T(x)=Q^T(x)*G^T(x)\pmod{x^{n-m+1}}$ ，消去了余数。

$F^T(x)=Q^T(x)*G^T(x)\pmod{x^{n-m+1}}$

变形得 $Q^T(x)=F^T(x)*G^T(x)^{-1}\pmod{x^{n-m+1}}$

这就可以直接用多项式求逆做了。求出 $Q^T(x)$ 后,再翻一次就得到 $Q(x)$ 。

$F(x)=Q(x)*G(x)+R(x)$ 变形得 $R(x)=F(x)-Q(x)*G(x)$ ，容易算出余数。

Code:

void mof(int *f,int *g,int n,int m){
  static int q[Maxn<<1],t[Maxn<<1];
  int L=n-m+1;
  rev(g,m);cpy(q,g,L);rev(g,m);
  rev(f,n);cpy(t,f,L);rev(f,n);
  invp(q,L);times(q,t,L,L);rev(q,L);
  times(g,q,n,n);
  for (int i=0;i<m-1;i++)
    g[i]=(f[i]-g[i]+mod)%mod;
  clr(g+m-1,L);
  cpy(f,q,L);clr(f+L,n-L);
}//F<-F/G;   G<-F%G.
int n,m,F[Maxn<<1],G[Maxn<<1];
int main()
{
  n=read()+1;m=read()+1;
  for (int i=0;i<n;i++)F[i]=read();
  for (int i=0;i<m;i++)G[i]=read();
  mof(F,G,n,m);
  print(F,n-m+1);print(G,m-1);
  return 0;
}

评测记录

4.多项式 $\ln,\exp$

两者均由麦克劳林级数定义：

$$\ln (F(x))=-\sum\limits_{i=1}\dfrac{(1-F(x))^i}{i}$$

\exp F(x)=\sum\limits_{i=0}\dfrac{F(x)^i}{i!}

为啥非得要用不友好的麦克劳林级数？因为这样就能做到仅使用加法和乘法来定义运算了。

经典的性质在这个定义下仍然成立，如 $\exp,\ln$ 互为逆运算， $\exp(\ln F(x)+\ln G(x))=F(x)*G(x)$ 等。

P4725 【模板】多项式对数函数（多项式 ln）

当然，可以直接使用定义式来计算答案，但是复杂度会十分糟糕。

考虑利用更加简明的性质，如 $\ln'(x)=\dfrac{1}{x}$ 。

题目有 $\ln(A(x))=B(x)$ ，两边同时求导。注意链式法则。

$\dfrac{d}{dx}\ln(A(x))=\dfrac{d}{dx}B(x)$

$\dfrac{dA(x)}{dx}\dfrac{d}{dA(x)}\ln(A(x))=B'(x)$

$\dfrac{A'(x)}{A(x)}=B'(x)$

同时积分得到 $B(x)=∫\!\left(A'(x)/A(x)\right)$

也就是说，一个求导，一个求逆，一个乘法，一个积分就好了。

Code:

void lnp(int *f,int m)
{
  static int g[Maxn<<1];
  cpy(g,f,m);
  invp(g,m);dao(f,m);
  times(f,g,m,m);
  jifen(f,m-1);
  clr(g,m);
}

评测记录

P4726 【模板】多项式指数函数（多项式 exp）

递推

有 $e^{F(x)}=G(x)$ ,则求导可得 $G'(x)=F'(x)G(x)$ (注意链式法则)

提取系数可得 $G'[n]=\sum\limits_{i=0}^nF'[i]G[n-i]$ 。

即 $(n+1)G[n+1]=\sum\limits_{i=0}^n(i+1)F[i+1]G[n-i]$ 。

每次可以 $O(n)$ 递推得到末项。

其实，由 $G(x)$ 也可推出 $F(x)$ ，即 $\ln$ 的递推算法。

由 $G'[n]=\sum\limits_{i=0}^nF'[i]G[n-i]$

$G'[n]=F'[n]+\sum\limits_{i=0}^{n-1}F'[i]G[n-i]$

$F'[n]=G'[n]-\sum\limits_{i=0}^{n-1}F'[i]G[n-i]$

$(n+1)F[n+1]=(n+1)G[n+1]-\sum\limits_{i=0}^{n-1}(i+1)F[i+1]G[n-i]$
倍增

和 $\ln$ 是逆运算 , 牛顿迭代大力搞。

回忆式子 $F(x)=F_*(x)-\dfrac{G(F_*(x))}{G'(F_*(x))}$

根据 $e^{A(x)}=B(x)$ ,设 $G(B(x))=\ln B(x)-A(x)$ ,那么这就是个牛顿迭代问题了。

求导得 $G'(B(x))=B(x)^{-1}$

$B(x)=B_*(x)-\dfrac{G(B_*(x))}{G'(B_*(x))}$

代入 $G$ 得到 $B(x)=B_*(x)-\dfrac{\ln B_*(x)-A(x)}{B_*(x)^{-1}}$

$B(x)=B_*(x)-(\ln B_*(x)-A(x))B_*(x)$

$B(x)=(1-\ln B_*(x)+A(x))B_*(x)$

根据上式倍增即可。

注意，必须保证原多项式 $A[0]=0$ ，此时 $B[0]=1$ 。否则在定义式中会出现无穷求和。

Code:

void exp(int *f,int m)
{
  static int s[Maxn<<1],s2[Maxn<<1];
  int n=1;for(;n<m;n<<=1);
  s2[0]=1;
  for (int len=2;len<=n;len<<=1){
    cpy(s,s2,len>>1);lnp(s,len);
  	for (int i=0;i<len;i++)
  	  s[i]=(f[i]-s[i]+mod)%mod;
  	s[0]=(s[0]+1)%mod;
  	times(s2,s,len,len);
  }cpy(f,s2,m);clr(s,n);clr(s2,n);
}

评测记录

分治FFT 做 $\exp$

观察上面的递推式，不难发现是个半在线卷积，分治FFT即可。

半在线卷积具体怎么做可见 : 半在线卷积小记

虽然复杂度是 $O(n\log^2 n)$ 的，但是由于常数较小，比大多数牛顿迭代快。代码和式子短，好写好调。

6.多项式快速幂

P5245 【模板】多项式快速幂

这里是 $O(n\log n)$ 的快速幂哦。

公式十分简单： $\large{A^k(x)=e^{\ln(A(x))*k}}$

( 也就是把乘法转化成指数上的加法 )

( 常数高达一个 $\log$ ,某些时候可以被 $O(n\log^2 n)$ 的暴力NTT代替 )

Code:

可以当做全套模板来用。

#include<algorithm>
#include<cstring>
#include<cstdio>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
const int _G=3,mod=998244353,Maxn=135000;
int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
inline void print(int *f,int len){
  for (int i=0;i<len;i++)
    printf("%d ",f[i]);
  puts("");
}
ll powM(ll a,ll t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  static ull f[Maxn<<1],w[Maxn]={1};
  tpre(n);
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }
    if (l==(1<<10))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
int _g1[Maxn<<1];
#define sav _g1
void times(int *f,int *g,int len,int lim)
{
  int n=1;while(n<len+len)n<<=1;
  cpy(sav,g,n);clr(sav+len,n-len);
  NTT(f,1,n);NTT(sav,1,n);
  px(f,sav,n);NTT(f,0,n);
  clr(f+lim,n-lim);clr(sav,n);
}
void invp(int *f,int m)
{
  int n;for (n=1;n<m;n<<=1);
  static int w[Maxn<<1],r[Maxn<<1];
  w[0]=powM(f[0]);
  for (int len=2;len<=n;len<<=1){
  	for (int i=0;i<(len>>1);i++)r[i]=w[i];
    cpy(sav,f,len);NTT(sav,1,len);
    NTT(r,1,len);px(r,sav,len);
    NTT(r,0,len);clr(r,len>>1);
    cpy(sav,w,len);NTT(sav,1,len);
  	NTT(r,1,len);px(r,sav,len);
  	NTT(r,0,len);
  	for (int i=len>>1;i<len;i++)
      w[i]=(w[i]*2ll-r[i]+mod)%mod;
  }cpy(f,w,m);clr(sav,n);clr(w,n);clr(r,n);
}
void dao(int *f,int m){
  for (int i=1;i<m;i++)
    f[i-1]=1ll*f[i]*i%mod;
  f[m-1]=0;
}
int inv[Maxn];
void Init(int lim){
  inv[1]=1;
  for (int i=2;i<=lim;i++)
    inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
void jifen(int *f,int m){
  for (int i=m;i;i--)
    f[i]=1ll*f[i-1]*inv[i]%mod;
  f[0]=0;
}
void lnp(int *f,int m)
{
  static int g[Maxn<<1];
  cpy(g,f,m);
  invp(g,m);dao(f,m);
  times(f,g,m,m);
  jifen(f,m-1);
  clr(g,m);
}
void exp(int *f,int m)
{
  static int s[Maxn<<1],s2[Maxn<<1];
  int n=1;for(;n<m;n<<=1);
  s2[0]=1;
  for (int len=2;len<=n;len<<=1){
    cpy(s,s2,len>>1);lnp(s,len);
  	for (int i=0;i<len;i++)
  	  s[i]=(f[i]-s[i]+mod)%mod;
  	s[0]=(s[0]+1)%mod;
  	times(s2,s,len,len);
  }cpy(f,s2,m);clr(s,n);clr(s2,n);
}
char str[Maxn];
int k,f[Maxn<<1],m;
int main()
{
  scanf("%d%s",&m,str);Init(m);
  for (int i=0;'0'<=str[i]&&str[i]<='9';i++)
    k=(10ll*k+str[i]-'0')%mod;
  for (int i=0;i<m;i++)f[i]=read();
  lnp(f,m);
  for(int i=0;i<m;i++)f[i]=1ll*f[i]*k%mod;
  exp(f,m);print(f,m);
  return 0;
}

评测记录

递推

设原多项式 $F(x)$ 的次数为 $k$ ,欲求 $F(x)^p$ 的前 $n$ 项。

对 $F(x)^{p+1}$ 求导,能得到 $(F(x)^{p+1})'=(p+1)F(x)^pF'(x)$

将 $F(x)^{p+1}$ 视作 $F(x)^p*F(x)$ 来求导，可求得 $(F(x)^{p+1})'=(F(x)^p)'F(x)+F(x)^pF'(x)$

拿这两个式子搞搞能得到 $p*F(x)^pF'(x)=(F(x)^p)'F(x)$

提取系数可得 : $p\sum\limits_{n-i\leq k-1}^{n}F'[n-i]F^p[i]=\sum\limits_{n-i\leq k}^{n}(i+1)F^p[i+1]F[n-i]$

挖出次数最高的项 :

$$p\sum\limits_{i=n-k+1}^{n}F'[n-i]F^p[i]-\sum\limits_{i=n-k}^{n-1}(i+1)F^p[i+1]F[n-i]=(n+1)F^p[n+1]F[0]$$

把 $i$ 的意义稍作修改,并整理。

$$F^p[n+1]=\dfrac{1}{(n+1)F[0]}\bigg(p\sum\limits_{i=0}^{k-1}(i+1)F[i+1]F^p[n-i]-\sum\limits_{i=1}^{k}F[i](n-i+1)F^p[n-i+1]\bigg)$$$$F^p[n+1]=\dfrac{1}{(n+1)F[0]}\bigg(p\sum\limits_{i=1}^{k}iF[i]F^p[n-i+1]-\sum\limits_{i=1}^{k}F[i](n-i+1)F^p[n-i+1]\bigg)$$$$F^p[n+1]=\dfrac{1}{(n+1)F[0]}\sum\limits_{i=1}^{k}(pi-n+i-1)F[i]F^p[n-i+1]$$

前面的和式枚举复杂度都是 $O(k)$ 的,一步步递推的复杂度就是 $O(nk)$ 的。

//默认f[0]=1
void powp(int *f,int k,int n,int p)
{
  static int g[Maxn]={1};
  for (int m=0;m<n-1;m++){
    for (int i=1;i<k;i++)
      g[m+1]=g[m+1]+(1ll*p*i-m+i-1)%mod*f[i]*g[m-i+1];
    g[m+1]=(ll)(g[m+1]+mod)*inv[m+1]%mod;
  }cpy(f,g,n);
}

这里也要求 $F[0]≠0$ ，否则需要做个次数平移。

( 令 $p=-1$ 就能得到求逆,令 $p=(mod-1)/2$ 能得到开根 )

注意，当 $k$ 较小的时候会比 Ln+EXP 快，后者只能跑 $10^5$ 的量级。

P5273 【模板】多项式幂函数 (加强版)

注意到其不保证满足 $A[0]=1$ ,无法直接套用前面的做法。

考虑把 $A(x)$ 写成 $B(x)*cx^p$ 的形式,让 $\ p=A(x)$ 末尾的空项个数 , $\ c=A(x)$ 最低非空项系数。

不难发现 $B[0]=1$ ,则 $B^k(x)$ 可求。

$A^k(x)=\big(B(x)*cx^k\big)^k=B^k(x)*c^kx^{pk}$

注意,常数 $c$ 的幂次要根据费马小定理,对 $mod-1$ 而非 $mod$ 取模。

求出 $B^k(x)$ 之后,乘一个单项式是容易的。

int main()
{
  scanf("%d%s",&m,str);
  for (int i=0;i<m;i++)f[i]=read();
  int p=0,c;while(!f[p])p++;c=powM(f[p]);
  for (int i=0;'0'<=str[i]&&str[i]<='9';i++){
    k=(10ll*k+str[i]-'0')%mod;
    k2=(10ll*k2+str[i]-'0')%(mod-1);
    if (1ll*k*p>m){
      for (int i=0;i<m;i++)printf("0 ");
      return 0;
    }
  }Init(m=m-p*k);
  for (int i=0;i<m;i++)f[i]=1ll*f[i+p]*c%mod;
  clr(f+m,p*k);lnp(f,m);
  for(int i=0;i<m;i++)f[i]=1ll*f[i]*k%mod;
  exp(f,m);
  for (int i=0;i<p*k;i++)printf("0 ");
  c=powM(c,mod-1-k2);
  for (int i=0;i<m;i++)f[i]=1ll*f[i]*c%mod;
  print(f,m);
  return 0;
}

评测记录

7.总结

讲了这么多，我们把几个算法的核心式子都集合一下。

NTT

$\large{F(ω^k_n)=Fl(ω^k_{n/2})+ω^k_nFR(ω^k_{n/2})}$

$\large{F(ω^{k+n/2}_n)=FL(ω^{k}_{n/2})-ω^k_nFR(ω^{k}_{n/2})}$

(其实主要是背代码)

可以通过ull以及预处理单位根优化。
多项式求逆

$G[0]$ 是 $F[0]$ 的逆元。

倍增公式 : $\large{G(x)=G^*(x)(2-F(x)G^*(x))}$

(带*的表示 $(mod\ x^{n/2})$ 的结果)
多项式牛顿迭代

已知 $G(x)$ 和 $G(F(x))=0$ ,求 $F(x)$ 。

倍增公式 : $\large{F(x)=F^*(x)-\frac{G(F^*(x))}{G'(F^*(x))}}(mod\ x^n)$
多项式开方

$B[0]$ 为 $A[0]$ 的二次剩余(同余系下的开方)

倍增公式 : $B(x)=\dfrac{A(x)+B^*(x)^2}{2B^*(x)}(mod\ x^n)$
多项式带余除法

定义一个 $n$ 次多项式翻转操作 $F^T(x)=x^nF(x^{-1})$ 。

$F^T(x)$ 的系数是 $F(x)$ 的反转。

把 $F(x)=Q(x)*G(x)+R(x)$ 里面的 $x$ 都换成 $x^{-1}$ ,等式两边同乘 $x^n$

得 $x^nF(x^{-1})=x^nQ(x^{-1})*G(x^{-1})+x^nR(x^{-1})$

变化为反转的形式，得 $F^T(x)=Q^T(x)*G^T(x)+x^{n-m+1}R^T(x)$

$\bmod\ x^{n-m+1}$ ，进而得到 $F^T(x)=Q^T(x)*G^T(x)(mod\ x^{n-m+1})$ ，消去了余数。

求出商之后余数易求。
多项式Ln

$(ln(x))'=\dfrac{1}{x}$

$ln'(A(x))A'(x)=B'(x)$

$B'(x)=\dfrac{A'(x)}{A(x)}=A'(x)A(x)^{-1}$

$\large{∫\!\left(A'(x)A(x)^{-1}\right)=B(x)}$

要求常数项为 $1$ ,则 $B[0]=0$ 。
多项式Exp

倍增公式 $\large{B(x)=(1-lnB^*(x)+A(x))B^*(x)}$ 。

要求常数项为 $0$ ，则 $B[0]=1$ 。
多项式快速幂

$\large{A^k(x)=e^{ln(A(x))*k}}$

附送一个 vector 板子 :

#include<algorithm>
#include<cstring>
#include<cstdio>
#include<vector>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
const int _G=3,mod=998244353,Maxn=1<<17|500;
using namespace std;
ll powM(ll a,ll t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn<<1],w[Maxn<<1];w[0]=1;
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }
    if (l==(1<<10))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
#define Poly vector<int>
Poly operator + (const Poly &A,const Poly &B)
{
  Poly C=A;C.resize(max(A.size(),B.size()));
  for (int i=0;i<B.size();i++)C[i]=(C[i]+B[i])%mod;
  return C;
}
Poly operator - (const Poly &A,const Poly &B)
{
  Poly C=A;C.resize(max(A.size(),B.size()));
  for (int i=0;i<B.size();i++)C[i]=(C[i]+mod-B[i])%mod;
  return C;
}
Poly operator * (const int c,const Poly &A)
{
  Poly C;C.resize(A.size());
  for (int i=0;i<A.size();i++)C[i]=1ll*c*A[i]%mod;
  return C;
}
int lim;
Poly operator * (const Poly &A,const Poly &B)
{
  static int a[Maxn<<1],b[Maxn<<1];
  for (int i=0;i<A.size();i++)a[i]=A[i];
  for (int i=0;i<A.size();i++)a[i]=A[i];
  cpy(a,&A[0],A.size());
  cpy(b,&B[0],B.size());
  Poly C;C.resize(min(lim,(int)(A.size()+B.size()-1)));
  int n=1;for(n;n<A.size()+B.size()-1;n<<=1);
  NTT(a,1,n);NTT(b,1,n);
  px(a,b,n);NTT(a,0,n);
  cpy(&C[0],a,C.size());
  clr(a,n);clr(b,n);
  return C;
}
void pinv(const Poly &A,Poly &B,int n)
{
  if (n==1)B.push_back(powM(A[0]));
  else if (n&1){
    pinv(A,B,--n);
    int sav=0;
    for (int i=0;i<n;i++)sav=(sav+1ll*B[i]*A[n-i])%mod;
    B.push_back(1ll*sav*powM(mod-A[0])%mod);
  }else {
    pinv(A,B,n/2);
    Poly sA;sA.resize(n);
    cpy(&sA[0],&A[0],n);
    B=2*B-B*B*sA;
    B.resize(n);
  }
}
Poly pinv(const Poly &A)
{
  Poly C;pinv(A,C,A.size());
  return C;
}
int inv[Maxn];
void Init()
{
  inv[1]=1;
  for (int i=2;i<=lim;i++)
    inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
Poly dao(const Poly &A)
{
  Poly C=A;
  for (int i=1;i<C.size();i++)
    C[i-1]=1ll*C[i]*i%mod;
  C.pop_back();
  return C;
}
Poly ints(const Poly &A)
{
  Poly C=A;
  for (int i=C.size()-1;i;i--)
    C[i]=1ll*C[i-1]*inv[i]%mod;
  C[0]=0;
  return C;
}
Poly ln(const Poly &A)
{return ints(dao(A)*pinv(A));}
void pexp(const Poly &A,Poly &B,int n)
{
  if (n==1)B.push_back(1);
  else if (n&1){
    pexp(A,B,n-1);n-=2;
    int sav=0;
    for (int i=0;i<=n;i++)sav=(sav+1ll*(i+1)*A[i+1]%mod*B[n-i])%mod;
    B.push_back(1ll*sav*inv[n+1]%mod);
  }else {
    pexp(A,B,n/2);
    Poly lnB=B;
    lnB.resize(n);lnB=ln(lnB);
    for (int i=0;i<lnB.size();i++)
      lnB[i]=(mod+A[i]-lnB[i])%mod;
    lnB[0]++;
    B=B*lnB;
    B.resize(n);
  }
}
Poly pexp(const Poly &A)
{
  Poly C;pexp(A,C,A.size());
  return C;
}
Poly F;
int main()
{
	int n;scanf("%d",&n);
  lim=n;Init();
	F.resize(n);
	for (int i=0;i<F.size();i++)scanf("%d",&F[i]);
	F=pexp(F);
	for (int i=0;i<F.size();i++)printf("%d ",F[i]);
	return 0;
}

8.高维多项式 / 多元幂级数

本部分未完。

先来介绍二元幂级数，即形如 $F(x,y)=\sum\limits_{i=0}\sum\limits_{j=0}F[i][j]x^iy^j$ 的含有两个元的幂级数。其系数是一个矩阵。

传统卷积

二元幂级数的乘法仍然定义为在各个维度上的加法卷积。

有 $(F*G)[n][m]=\sum\limits_{i=0}^n\sum\limits_{j=0}^mF[i][j]*G[n-i][m-j]$。

二元幂级数的“系数”可以仍然是幂级数，如 $(G[n])(y)=[x^n]G(x,y)=\sum\limits_{i=0}G[n][i]y^i$

也有 $\big((F*G)[n]\big)(y)=\sum\limits_{i=0}^nF[i](y)*G[n-i](y)$

可以先对 $y$ 的维度做 $\rm DFT$ ，得到点值。

然后上式中 $F[i](y)*G[n-i](y)$ 就变成了普通的“数值”乘法，整个式子也就变成了 $x$ 上的普通加法卷积。

这可以得到二元幂级数卷积的过程 :

先对行做 $\rm DFT$ ，再对列做 $\rm DFT$ ，得到点值，点乘之后先对列做 $\rm IDFT$ ，再对行做 $\rm IDFT$ 。

也可以通过二元函数插值来解释。

全家桶

有了卷积，我们再来推导多项式全家桶。

求逆仍然可以使用倍增公式 $G(x,y)=G^*(x,y)(2-F(x,y)G^*(x,y))$ ，因为平方法的普适性很强。

$\ln,\exp$ 仍然可以使用级数定义。

咕。

维度爆炸

$\text{stO EI Orz}$ ,下面的内容大多都是复读。

设多项式共有 $k$ 个维度，界分别为 $n_1,n_2,...,n_k$ 。设 $N=\prod_{i=1}^kn_i$ 。

在高维卷积中，由于每一维度上次数都会翻倍，所以实际计算系数总量是 $O(N2^k)$ 级别的。当 $k$ 较大时，复杂度不尽如人意。

下面给出一个 $O(kN\log N)$ 计算 $k$ 维多项式乘法的方法。鉴于 $n_i$ 一般至少为 $2$ ，这里的 $k$ 应最多是 $O(\log N)$ 级别。

将一个下标 $(i_1,i_2,...,i_k)$ 编码为一个 $(n_1,n_2,...,n_k)$ 进制数。

即令 $i=\sum\limits_{t=1}^ki_t\prod_{r=1}^{t-1}n_r=i_1+i_2*n_1+i_3*n_1*n_2+...+i_k*n_1*n_2...*n_{k-1}$。

这样，就把多维映射到了一维上。

此时，下标 $(i_1,i_2,...,i_k)$ 的加法正对应着大数 $i$ 的加法，但要将超出范围(产生进位)的贡献去除。

考虑设置一个合适的占位多项式来辅助判定进位，分流贡献。考虑占位函数 $χ(n)$ ，满足 $i+j$ 不进位当且仅当 $χ(i)+χ(j)=χ(i+j)$ 。

这样，我们计算二元幂级数 $\sum\limits_{i=0}F[i]x^iy^{χ(i)}$ 的卷积，然后利用第二维去除无用贡献即可。

现在剩下的问题就是构造出一个合适的 $χ$ 函数。

模仿子集卷积，不难想到 $χ(i)=\sum\limits_{t=1}^ki_t$ ，但是这个函数的值域太大，会导致复杂度退化。

然后就是 $\rm EI$ 大神古今一绝的构造了。

注意到，令 $P=\prod_{r=1}^{t-1}n_r$ ， $i+j$ 在第 $t$ 位进位当且仅当 $\left\lfloor \dfrac{i}{P}\right\rfloor+\left\lfloor \dfrac{j}{P}\right\rfloor+1=\left\lfloor \dfrac{i+j}{P}\right\rfloor$

令 $χ(i)=\sum\limits_{t=1}^{k-1}\left\lfloor \dfrac{i}{\prod_{r=1}^{t}n_r}\right\rfloor=\left\lfloor \dfrac{i}{n_1}\right\rfloor+\left\lfloor \dfrac{i}{n_1n_2}\right\rfloor+...+\left\lfloor \dfrac{i}{n_1n_2...n_k}\right\rfloor$

这样，虽然单个 $χ(i)$ 的值仍然很大，但是由于 $\left\lfloor\dfrac{i+j}{P}\right\rfloor-\Bigg(\left\lfloor\dfrac{i}{P}\right\rfloor+\left\lfloor\dfrac{j}{P}\right\rfloor\Bigg)\in\{0,1\}$ ，所以有 $χ(i+j)-\Big(χ(i)+χ(j)\Big)\in[0,k)$ 。

这样，只需记录模 $y^k-1$ 的循环卷积即可得到所有信息。

实现时，第二维可以单独封装成一个结构体。这样也可以直接套用一元牛顿迭代。

板子：[数学记录]Uoj#596. 【集训队互测2021】三维立体混元劲

9.任意模数多项式乘法

P4245 【模板】任意模数NTT

有的时候题目给的膜数并不是 $998244353$ 等和蔼的数字，而可能是 $10^9+7$ 这种（~~看似平淡无奇实则暗藏杀机~~）。

那么，我们苦苦研究的 NTT 就都没有用武之地了吗？

任意模数多项式乘法原理 : 假设卷积前，每个数都小于 $mod$ 。

卷积后，能产生的最大的数也就是 $mod*mod*len$ ，实战应用中一般约为 $10^9*10^9*10^5=10^{23}$ 。

我们只要弄一个精度够高的多项式乘法，就能做到不丢精度。

一种想法是 : 跑九次 NTT (三次卷积)，把答案用中国剩余定理合并，精度可达 $10^{26}$ 。

这种做法常数非常大，而且 $10^{26}$ 用 long long 存不下，还需要一些黑科技辅助，所以不推荐。

另一种想法是 : 拆系数 FFT ，又称为 MTT 。

把一个数拆成 $a*2^{15}+b$ 的形式，那么 $a,b<=2^{15}$ 。

把 $a$ 和 $b$ 弄成两个多项式，这样相乘的值域是 $n*\sqrt{mod}*\sqrt{mod}≈10^{14}$

那么 $c_1*c_2=(a_1*2^{15}+b_1)(a_2*2^{15}+b_2)$

$=a_1a_2*2^{30}+(a_1b_2+a_2b_1)2^{15}+b_1b_2$

这样做需要 $4$ 次多项式乘法。12次 FFT? 常数爆炸!

冷静分析 : 每个多项式只用插值一次，共 $4$ 次。

点乘复杂度忽略，最后得到的四个多项式各插值一次，共 $4$ 次。

这样就是 $8$ 次 FFT ，常数还是爆炸……

我们考虑推推式子来优化:

我们有四个多项式 $A1,A2,B1,B2$ ,求这些多项式的两两乘积。

考虑到 $(a+bi)*(c+di)=a(c+di)*bi(c+di)=ac-bd+(ad+bc)i$

那么我们设复多项式 $P=A1+iA2;\ Q=B1+iB2;$

那么 $T1=P*Q=A1B1-A2B2+(A1B2+A2B1)i$

我们又设 $P'=A1-iA2$

那么 $T2=P'*Q=A1B1+A2B2+(A1B2-A2B1)i$

$T1+T2=2(A1B1+iA1B2)$ ,这样我们就求出了其中两个，减一减就能得到另外两个。

总的 FFT 次数是 : $3$ 次DFT $+2$ 次IDFT $=5$ 次.

不过，值域 $10^{14}$ 有点吃紧，long double信仰跑吧。

提高精度的方法：预处理单位根。这样可以使得每个单位根都仅用 $O(\log)$ 次乘法算出。具体方法见代码。

封装版 :

#include<algorithm>
#include<cstdio>
#include<cmath>
#define ld /*long*/ double
#define ll long long
#define Maxn 135000
const ld Pi=acos(-1);
using namespace std;
inline int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
struct CP{
  ld x,y;
  CP operator + (const CP& B) const
  {return (CP){x+B.x,y+B.y};}
  CP operator - (const CP& B) const
  {return (CP){x-B.x,y-B.y};}
  CP operator * (const CP& B) const
  {return (CP){x*B.x-y*B.y,x*B.y+y*B.x};}
};
int mod,tr[Maxn<<1];
void FFT(CP *f,int op,int n)
{
  static CP w[Maxn]={(CP){1.0,0.0}};
  for (int i=0;i<n;i++)
    if (i<tr[i])swap(f[i],f[tr[i]]);
  for(int l=1;l<n;l<<=1){
    CP tG=(CP){cos(Pi/l),sin(Pi/l)*op};
    for (int i=l-2;i>=0;i-=2)w[i+1]=(w[i]=w[i>>1])*tG;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        CP sav=w[p]*f[k|l|p];
        f[k|l|p]=f[k|p]-sav;
        f[k|p]=f[k|p]+sav;
      }
  }
}
void times(int *f,int *g,int n,int lim)
{
  static CP P1[Maxn<<1],P2[Maxn<<1],Q[Maxn<<1];
  for (int i=0,sav;i<n;i++){
    P1[i]=(CP){f[i]>>15,f[i]&32767};
    P2[i]=(CP){f[i]>>15,-(f[i]&32767)};
    Q[i]=(CP){g[i]>>15,g[i]&32767};
  }
  for (int i=1;i<n;i++)
    tr[i]=tr[i>>1]>>1|((i&1)?n>>1:0);
  FFT(P1,1,n);FFT(P2,1,n);FFT(Q,1,n);
  for (int i=0;i<n;i++){
    Q[i].x/=n;Q[i].y/=n;
    P1[i]=P1[i]*Q[i];
    P2[i]=P2[i]*Q[i];
  }FFT(P1,-1,n);FFT(P2,-1,n);
  for (int i=0;i<lim;i++){
    ll a1b1=0,a1b2=0,a2b1=0,a2b2;
    a1b1=(ll)floor((P1[i].x+P2[i].x)/2+0.49)%mod;
    a1b2=(ll)floor((P1[i].y+P2[i].y)/2+0.49)%mod;
    a2b1=((ll)floor(P1[i].y+0.49)-a1b2)%mod;
    a2b2=((ll)floor(P2[i].x+0.49)-a1b1)%mod;
    f[i]=((((a1b1<<15)+(a1b2+a2b1))<<15)+a2b2)%mod;
    if (f[i]<0)f[i]+=mod;
  }
}
int n,m,f[Maxn<<1],g[Maxn<<1];
int main()
{
  scanf("%d%d%d",&n,&m,&mod);n++;m++;
  for (int i=0;i<n;i++)f[i]=read()%mod;
  for (int i=0;i<m;i++)g[i]=read()%mod;
  int len=1;
  for (m=m+n-1;len<m;len<<=1);
  times(f,g,len,m);
  for (int i=0;i<m;i++)printf("%d ",f[i]);
  return 0;
}

10.后记

在寒假的最后几天,我一边补着如山的寒假作业,一遍写着这篇文章。

如果有纰漏的话请私信指正,如果有问题的话也欢迎提问。

欲知后事如何,请看下节 : 多项式计数杂谈

NTT与多项式全家桶

-1.原根的引入

0.NTT

$\frac{1}{2}$ .注意事项

1.多项式四则运算

2.多项式牛顿迭代(及求导/积分/复合)

3.多项式带余除法

4.多项式 $\ln,\exp$

6.多项式快速幂

7.总结

NTT

多项式求逆

多项式牛顿迭代

多项式开方

多项式带余除法

多项式Ln

多项式Exp

多项式快速幂

8.高维多项式 / 多元幂级数

9.任意模数多项式乘法

10.后记

状态

开发

支持

关于

NTT与多项式全家桶

-1.原根的引入

0.NTT

12\frac{1}{2}21​.注意事项

1.多项式四则运算

2.多项式牛顿迭代(及求导/积分/复合)

3.多项式带余除法

4.多项式 ln⁡,exp⁡\ln,\expln,exp

6.多项式快速幂

7.总结

NTT

多项式求逆

多项式牛顿迭代

多项式开方

多项式带余除法

多项式Ln

多项式Exp

多项式快速幂

8.高维多项式 / 多元幂级数

9.任意模数多项式乘法

10.后记

状态

开发

支持

关于

还没有账户？

登录

$\frac{1}{2}$ .注意事项

4.多项式 $\ln,\exp$