NTT与多项式全家桶

command_block

2019-01-31 13:47:34

Personal

(\scriptsize\text{由于文章前后内容有所差别,评论区已于2019.10.28清屏})

(\scriptsize\text{2020.03.27 : 发觉了之前的naive,开始分拆封装卡常数})

(\scriptsize\text{2020.05.11 : 更新了更快的NTT})

(\scriptsize\text{2021.01.10 : 修改了若干错误,加入了更严谨的描述,微调顺序;更新高维多项式})

(\scriptsize\text{2021.01.11 : 更快的求逆})

(\scriptsize\text{2021.02.24 : vector板子})

书接上回 : 傅里叶变换(FFT)学习笔记

学完了FFT乘法之后这个坑就鸽了半年……

这里的全家桶指 {求逆+牛顿迭代+带余除法+Ln+Exp+快速幂+MTT+高维多项式}

此外,插值求值请见 从拉插到快速插值求值

如果对生成函数运算的定义方式与合法性感兴趣,可见 : rqy : 浅谈 OI 中常用的一些生成函数运算的合法与正确性

-1.原根的引入

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

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

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

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

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

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

我们来看看原根的定义。

众所周知,对于素数 p ,其模意义同余系 F_p 中恰有 1,2,3,...p-1p-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^rs 为奇数。

随机出一个二次非剩余 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 变换的时间为 FFT75%~80%.

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

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

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

一些额外的技巧:

在不同的题目中,可以加速 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) 变换的线性性可以减少运算。

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

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.多项式四则运算

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

多项式的加减定义为 :

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]

使用 NTTO(n\log n) 内完成计算。

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

接下来就是多项式乘法逆元,定义为 :

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

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

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

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) 即可得到上式,利用整体法不难理解。

用于解决下列问题:

已知函数 GG(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}} 级别。

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^?} 来消去某些项。很明显,这种操作只能 消去高次的项,而留下低次的项

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

首先有 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)

**Code:** ```cpp 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; } ``` [评测记录](https://www.luogu.com.cn/record/44844107) # 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)

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即可。

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

虽然复杂度是 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 的量级。

注意到其不保证满足 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.总结

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

附送一个 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)) ,因为平方法的普适性很强。

咕。 - **维度爆炸** [$\text{stO EI Orz}$](https://www.luogu.com.cn/blog/EntropyIncreaser/hello-multivariate-multiplication) ,下面的内容大多都是复读。 设多项式共有 $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}

ab 弄成两个多项式,这样相乘的值域是 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

总的 FFT 次数是 : $3$ 次DFT $+2$ 次IDFT $=5$ 次. 不过,值域 $10^{14}$ 有点吃紧,`long double`信仰跑吧。 提高精度的方法 : 预处理单位根。这样可以使得每个单位根都仅用 $O(\log)$ 次乘法算出。具体方法见代码。 封装版 : ```cpp #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.后记 在寒假的最后几天,我一边补着如山的寒假作业,一遍写着这篇文章。 如果有纰漏的话请私信指正,如果有问题的话也欢迎提问。 欲知后事如何,请看下节 : [多项式计数杂谈](https://www.luogu.com.cn/blog/command-block/sheng-cheng-han-shuo-za-tan)