NTT与多项式全家桶
command_block
2019-01-31 13:47:34
($\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)学习笔记](https://www.luogu.org/blog/command-block/fft-xue-xi-bi-ji)
学完了FFT乘法之后这个坑就鸽了半年……
这里的**全家桶**指 {求逆+牛顿迭代+带余除法+Ln+Exp+快速幂+MTT+高维多项式}
此外，插值求值请见 [从拉插到快速插值求值](https://www.luogu.org/blog/command-block/zong-la-cha-dao-kuai-su-cha-zhi-qiu-zhi)
如果对生成函数运算的定义方式与合法性感兴趣，可见 ： [rqy : 浅谈 OI 中常用的一些生成函数运算的合法与正确性](https://rqy.moe/Math/gf_correct/)
# -1.原根的引入
FFT 对单位根的依赖，决定该算法必须采用浮点数，进而引发精度问题。
糟糕的是，数学家们已经证明，在复数域 $\mathbb C$ 中，单位根是唯一一类满足要求的数。
大部分的计数题都是在模意义下完成的，我们自然希望为 FFT 找一个模意义下的替代品。
考虑引入模意义同余系 $F_p$ 中单位根的类似物 : **原根**。
可能有帮助的文章 : [原根&离散对数相关](https://www.luogu.com.cn/blog/command-block/yuan-gen-li-san-dui-shuo-xiang-guan)
我们先罗列一下 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$ 即为一个很有代表性的模数。
熟练的选手会记下一些常用模数的原根。
懒人福利 : 安利[大佬的原根表](http://blog.miskcoo.com/2014/07/fft-prime-table)。
至于如何寻找原根，可见上面拉链的那篇文章。
但是，一定非要原根不可吗?
摆脱原根的方法来自 [Uoj : hly1204](https://hly1204.blog.uoj.ac/blog/6422)
假如我们已经知道 $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评测记录](https://www.luogu.org/record/25876484)
```cpp
#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出现!
还是原来的配方,还是熟悉的[味道](https://www.luogu.org/problemnew/show/P3803)。
**Code**
```cpp
#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;
}
```
[评测记录](https://www.luogu.com.cn/record/32183302)
未经毒瘤卡常时,单个 `NTT` 变换的时间为 `FFT` 的 `75%~80%`.
其中这一句话改变得很大:
```cpp
ll tG=powM(op==1 ? G:invG,(mod-1)/p);
```
即 $\omega_p^1=g^{(mod-1)/p}$ ,再也不用三角函数啦!
一些额外的技巧:
- 预处理单位根，这样可以省去部分乘法操作，并访问一段连续内存。
- 注意到只会有加减法操作，可以使用 `ull` 存储，能耐受 `18*mod*mod` 的范围，在常规范围下可以不取模。
模板题范围较大,在中间取一次模即可。
在不同的题目中，可以加速 `15%~40%` 不等。
```cpp
#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** (仅模板)
```cpp
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("");
}
```
- **多项式基本操作**
```cpp
#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 中常用的一些生成函数运算的合法与正确性](https://rqy.moe/Math/gf_correct/)。
简化地讲，我们将以**多项式加法**和**加法卷积**来定义所有运算。
- **界**
在大多数题目中，我们只对多项式的前若干项感兴趣，此时我们为所有运算设定一个次数上界，即 $\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 【模板】多项式乘法逆](https://www.luogu.org/problemnew/show/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** : (配合前面的封装食用)
```cpp
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`。
检查的方法 : 随机一个多项式求逆两次看看是不是自己。
[评测记录](https://www.luogu.com.cn/record/32183046)
- 卡常后的倍增
众所周知 $\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$ 项，正好是我们不需要的。
```cpp
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);
}
```
[评测记录](https://www.luogu.com.cn/record/44843813) (可以加速两倍多)
# 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:**
```cpp
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$ 一般较简单，而且是个常量(可以手动数学分析)。
证明可在 [多项式计数杂谈](https://www.luogu.com.cn/blog/command-block/sheng-cheng-han-shuo-za-tan) 中查看。
**结论** ：$\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 【模板】多项式开根](https://www.luogu.org/problemnew/show/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:**
```cpp
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);
}
```
[评测记录](https://www.luogu.com.cn/record/44844033)。
# 3.多项式带余除法
[P4512 【模板】多项式除法](https://www.luogu.org/problemnew/show/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:**
```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）](https://www.luogu.org/problemnew/show/P4725)
当然，可以直接使用定义式来计算答案，但是复杂度会十分糟糕。
考虑利用更加简明的性质，如 $\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:**
```cpp
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);
}
```
[评测记录](https://www.luogu.org/recordnew/show/18245545)
[P4726 【模板】多项式指数函数（多项式 exp）](https://www.luogu.org/problemnew/show/P4726)
- **递推**
有 $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:**
```cpp
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);
}
```
[评测记录](https://www.luogu.com.cn/record/44844364)
- 分治FFT 做 $\exp$
观察上面的递推式，不难发现是个半在线卷积，分治`FFT`即可。
半在线卷积具体怎么做可见 : [半在线卷积小记](https://www.luogu.com.cn/blog/command-block/ban-zai-xian-juan-ji-xiao-ji)
虽然复杂度是 $O(n\log^2 n)$ 的，但是由于常数较小，比大多数牛顿迭代快。代码和式子短，好写好调。
# 6.多项式快速幂
[P5245 【模板】多项式快速幂](https://www.luogu.org/problemnew/show/P5245)
这里是 $O(n\log n)$ 的快速幂哦。
公式十分简单：$\large{A^k(x)=e^{\ln(A(x))*k}}$
( 也就是把乘法转化成指数上的加法 )
( 常数高达一个 $\log$ ,某些时候可以被 $O(n\log^2 n)$ 的暴力NTT代替 )
**Code:**
可以当做全套模板来用。
```cpp
#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;
}
```
[评测记录](https://www.luogu.com.cn/record/44844420)
- **递推**
设原多项式 $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)$的。
```cpp
//默认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 【模板】多项式幂函数 (加强版)](https://www.luogu.com.cn/problem/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)$ 之后,乘一个单项式是容易的。
```cpp
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;
}
```
[评测记录](https://www.luogu.com.cn/record/34015941)
# 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` 板子 :
```cpp
#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}$](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】三维立体混元劲](https://www.luogu.com.cn/blog/command-block/post-shuo-xue-ji-lu-uoj596-ji-xun-dui-hu-ce-2021-san-wei-li-ti-hun-y)
# 9.任意模数多项式乘法
[P4245 【模板】任意模数NTT](https://www.luogu.org/problemnew/show/P4245)
有的时候题目给的膜数并不是 $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)$ 次乘法算出。具体方法见代码。
封装版 :
```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)