多项式总结
Great_Influence
2018-05-11 18:37:51
此处多项式算法，主要针对基于快速傅里叶变换 $(FFT)$ 或者快速数论变换 $(NTT)$ 的多项式演算，主要用来解决各种母函数、排列组合、字符串题目。此处仅介绍相关底层算法。文章可能不定期更新。
[TOC]
## 1.多项式乘法
[点我](https://www.luogu.org/problemnew/show/P3803)
一切算法的基础。直接做为 $O(n^2)$ ，利用 $FFT$ 或者 $NTT$ 来实现可以做到 $O(n\log n)$。
配合预处理单位根可以提升精度和速度。
#### 代码:
```cpp
inline void FFT(comp*F,int tl)
{
rep(i,0,Len)if(i<rev[i])swap(F[rev[i]],F[i]);
for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
for(register int j=0;j<Len;j+=i)rep(k,0,ii)
{
comp x=F[j+k+ii]*g[t][k];
F[j+k+ii]=F[j+k]-x;
F[j+k]=F[j+k]+x;
}
if(tl==-1)
{
reverse(F+1,F+Len);
rep(i,0,Len)F[i]=F[i]/Len;
}
}
inline void NTT(int*F,int typ)
{
Rep(i,1,Len-1)if(i<rev[i])swap(F[i],F[rev[i]]);
for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
for(register int j=0;j<Len;j+=i)rep(k,0,ii)
{
register int tt=(uint64)F[j+k+ii]*g[t][k]%mod;
F[j+k+ii]=ad(F[j+k],mod-tt);
F[j+k]=ad(F[j+k],tt);
}
if(typ==-1)
{
reverse(F+1,F+Len);
register uint64 invn=power(Len,mod-2);
rep(i,0,Len)F[i]=invn*F[i]%mod;
}
}
int X[MAXN],Y[MAXN];
inline void mul(int *a,int *b,int *c,int lenl,int lenr)
{
if(lenl<=300/lenr)
{
Rep(i,0,lenl+lenr)c[i]=0;
Rep(i,0,lenl)Rep(j,0,lenr)c[i+j]=(c[i+j]+a[i]*b[j])%mod;
return;
}
for(Len=2;Len<=lenl+lenr;Len<<=1);
Rep(i,0,lenl)X[i]=a[i];
Rep(i,0,lenr)Y[i]=b[i];
rep(i,lenl+1,Len)X[i]=0;
rep(i,lenr+1,Len)Y[i]=0;
NTT(X,1),NTT(Y,1);
rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
NTT(X,-1);
Rep(i,0,lenl+lenr)c[i]=X[i];
}
```
注意 $FFT$ 的复数类最好自己手写，比系统自带的快多了。
以下几个操作大多基于取模来实现，存在一个特殊的模数，满足其减 $1$ 后存在较多 $2$ 的因数。无说明一般为 $998244353$ 。
## 2.多项式求逆
[点我](https://www.luogu.org/problemnew/show/P4238)
也算是基础吧。
#### 推导过程:
设要求 $F$ 关于 $x^n$ 的逆，则
设
$f_i(x)\equiv F^{-1}(x)\pmod{x^{2^i}}$
假设已知 $f_i$ ，要求 $f_{i+1}$ ，那么
$F(x)f_i(x)\equiv F(x)f_{i+1}(x)\pmod{x^{2^i}}$
$[f_i(x)-f_{i+1}(x)]\equiv 0\pmod{x^{2^i}}$
两边平方，则
$[f_i(x)-f_{i+1}(x)]^2\equiv 0\pmod{x^{2^{i+1}}}$
$f_i^2(x)-2f_i(x)f_{i+1}(x)+f_{i+1}^2(x)\equiv 0\pmod{x^{2^{i+1}}}$
同乘$F$,得
$f_{i+1}(x)-2f_i(x)+F(x)f_i^2(x)\equiv 0\pmod{x^{2^{i+1}}}$
$f_{i+1}(x)\equiv 2f_i(x)-F(x)f_i^2(x)\pmod{x^{2^{i+1}}}$
利用这个公式直接倍增就好了。时间复杂度利用主定理可以证明为 $O(n\log n)$ 。
#### 代码:
```cpp
inline void Inv(int*F,int*G,int ln)
{
Iv[0]=power(F[0],mod-2);
rep(i,0,ln)X[i]=Y[i]=0;
for(register int Ln=2;Ln<=ln;Ln<<=1)
{
rep(i,0,Ln)X[i]=F[i];
rep(i,0,(Ln>>1))Y[i]=Iv[i];
Len=Ln,calrev();
NTT(X,1),NTT(Y,1);
rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
NTT(X,-1);
rep(i,0,(Ln>>1))X[i]=0;
NTT(X,1);
rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
NTT(X,-1);
rep(i,(Ln>>1),Ln)Iv[i]=mod-X[i];
}
rep(i,0,ln)G[i]=Iv[i];
}
```
## 3.多项式除法
[点我](https://www.luogu.org/problemnew/show/P4512)
已知函数$F$和$G$，求$Q$和$R$，使得
$F=QG+R$
保证 $Deg(R)<Deg(G)<Deg(F)$ 。
#### 推导过程:
先考虑一种操作:
$F_R(x)=x^{Deg(F)}F(\frac{1}{x})$
拆开发现该操作的作用就是将$F$的系数反位。
然后开始化式子。
为了简便，设$Deg(F)=n,Deg(G)=m$，那么可知$Deg(Q)=n-m,Deg(R)=m-1$。
$\because F(x)=G(x)Q(x)+R(x)$
$\therefore F(\frac{1}{x})=G(\frac{1}{x})Q(\frac{1}{x})+R(\frac{1}{x})$
$\therefore x^nF(\frac{1}{x})=x^nG(\frac{1}{x})Q(\frac{1}{x})+x^nR(\frac{1}{x})$
$\therefore x^nF(\frac{1}{x})=x^mG(\frac{1}{x})x^{n-m}Q(\frac{1}{x})+x^{n-m+1}*x^{m-1}R(\frac{1}{x})$
$\therefore F_R(x)=G_R(x)Q_R(x)+x^{n-m+1}*R_R(x)$
将这个式子放在模$n-m+1$次方下，可以发现$R_R$的所有系数正好全部被模掉，失去影响。
$\therefore F_R(x)=G_R(x)Q_R(x)\pmod{x^{n-m+1}}$
$\therefore Q_R(x)=F_R(x)G_R^{-1}(x)\pmod{x^{n-m+1}}$
多项式求逆即可。至于$R$，直接$R=F-QG$就可以求出。时间复杂度 $O(n\log n)$ 。
#### 代码:
```cpp
static int X[MAXN],Y[MAXN];
inline void Div(int *a,int n,int *b,int m)
{
memcpy(X,a,sizeof X);memcpy(Y,b,sizeof Y);
reverse(b,b+m+1);reverse(X,X+n+1);
Rep(i,n-m+1,n)X[i]=0;
Inv(b,n-m+1);
while(Len<=(n<<2))Len<<=1;
calrev(Len);
mul(X,b);reverse(X,X+n-m+1);
Rep(i,n-m+1,n)X[i]=0;
mul(Y,X);
Rep(i,0,m-1)b[i]=ad(a[i],mod-Y[i]);
memcpy(a,X,sizeof X);
}
```
## 4.多项式开方
[点我](https://www.luogu.org/problemnew/show/CF438E)
开始有点麻烦了。依旧采用倍增来求。
已知$F$，求$G$使得$G^2\equiv F\pmod{x^n}$
假设求出 $T_0^2\equiv F\pmod{x^{t}}$ ，要求 $T^2\equiv F\pmod{x^{2t}}$ 。
有
$T_0^2-T^2\equiv 0\pmod{x^t}$
$(T_0+T)(T_0-T)\equiv 0\pmod{x^t}$
因为存在平方，所以存在2个解。设 $T-T_0\equiv 0\pmod{x^{t}}$ ，则再平方有
$T'^2-2T_T+T^2\equiv 0\pmod{x^{2t}}$
又因为 $T^2\equiv F\pmod{x^{2t}}$ ，则
$T'^2-2TT_0+F\equiv 0\pmod{x^{2t}}$
$T\equiv \frac{T_0^2+F}{2T_0}\pmod{x^{2t}}$
直接倍增计算即可。时间复杂度依旧 $O(n\log n)$ 。
#### 代码:
```cpp
inline void Sqrt(int*F,int*G,int ln)
{
G[0]=1;
for(register int Ln=2;Ln<=ln;Ln<<=1)
{
mul(G,G,ExX,Ln>>1,Ln>>1);
rep(i,0,Ln>>1)ExX[i]=ad(mod-ExX[i+(Ln>>1)],F[i+(Ln>>1)]);
rep(i,Ln>>1,Ln)ExX[i]=0;
Inv(G,ExY,Ln>>1);
Len=Ln,calrev();
NTT(ExX,1),NTT(ExY,1);
rep(i,0,Ln)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
NTT(ExX,-1);
rep(i,0,Ln>>1)G[i+(Ln>>1)]=(ll)ExX[i]*inv2%mod;
}
}
```
## 5.多项式积分和求导
~~这就不要例题了吧...~~
设$F(x)=\sum\limits_{i=0}^n a_ix^i$，则
定义:
$\int F(x)=\sum\limits_{i=0}^n \frac{a_i}{i+1}x^{i+1}$
$\frac{dF}{dx}=\sum\limits_{i=1}^n ia_ix^{i-1}$
直接$O(n)$计算就可以了。
#### 代码:
```cpp
inline void Direv(int *a,int n)
{Rep(i,1,n)a[i-1]=(ll)a[i]*i%mod;a[n]=0;}
inline void Inter(int *a,int n)
{Repe(i,n,1)a[i]=(ll)a[i-1]*power(i,mod-2)%mod;a[0]=0;}
```
## 6.多项式取 $\ln$
[模板题点我](https://www.luogu.org/problemnew/show/P4725)
[做完模板可以点我](https://www.lydsy.com/JudgeOnline/problem.php?id=3456)
$\ln$ 的定义为 $\exp$ 的逆变换。具体定义请参考下面的 $\exp$ 。
定义式为:
$$
\ln(1-x)=-\sum_{i=1}^\infty \frac{x^i}{i}
$$
存在推导式:
$$\ln F=\int \frac{F'}{F}\equiv \int F'F^{-1}$$
直接计算即可。
#### 代码:
```cpp
inline void Ln(int *a,int n)
{
memcpy(X,a,sizeof X);
Inv(X,n);Direv(a,n);
while(Len<=(n<<2))Len<<=1;
calrev(Len);
mul(a,X);
Inter(a,n);
Rep(i,n+1,Len)a[i]=0;
}
```
## 7.多项式牛顿迭代
题目请见 $WC2019$ 汪乐平的课件。
采用倍增实现。
已知$G$，求多项式$F$使得$G(F(x))= 0$。
设已知$F_0$满足$G(F_0(x))\equiv 0\pmod{x^t}$，要求$G(F(x))\equiv 0\pmod{x^{2t}}$。
存在递推式
$F(x)≡F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))} \pmod {x^{2t}}$
其中，$G'(F(x))=\frac{\mathbb{d}G}{\mathbb{d}F}$，即 $G$ 对 $F$ 的导数。
直接计算即可。
### 7.1推导过程:
函数 $G$ 对 $F$ 进行泰勒展开，得:
$$G(F)=G(F_0)+G'(F_0)(F-F_0)+G''(F_0)(F-F_0)^2+\cdots$$
因为$F\equiv F_0\pmod{x^t}$，因此$F-F_0\equiv 0\pmod{x^t}$，因此$(F-F_0)^k\equiv 0\pmod{x^{2t}},k\ge 2$。
因此后面几项在$\bmod x^{2t}$意义下全部为$0$，直接忽略。
因此得到关系式
$$G(F)=G(F_0)+G'(F_0)(F-F_0)$$
整理即得上式。
(代码依 $G$ 不同而不同)
### 7.2.其他式子利用牛顿迭代推导的方式
#### 7.2.1 求逆
考虑对方程$F(x)G(x)-1=0$ 作牛顿迭代，得:
$\delta(F)=GF-1$
$\delta'(F)=G$
可以得到递推式:
$F=F_0-\frac{\delta(F_0)}{\delta'(F_0)}$
$=F_0-\frac{G_0F_0-1}{G_0}$
$=F_0-F_0(G_0F_0-1)$
$=2F_0-G_0F_0^2$
#### 7.2.2 开方
考虑对方程$F^2-G=0$作牛顿迭代，得:
$\delta(F)=F^2-G$
$\delta'(F)=2F$
可以得到递推式:
$F=F_0-\frac{\delta(F_0)}{\delta'(F_0)}$
$=F_0-\frac{F_0^2-G_0}{2F_0}$
$=\frac{F_0^2+G_0}{2F_0}$
#### 7.2.3 $\exp$
见下文。
### 7.3 牛顿法的通用优化技巧
使用牛顿法，我们通常需要多次求出$F_0-\frac{G(F_0)}{G'(F_0)}$。有时为了方便，我们会将前面的$F_0$和后面部分共同计算。
可以发现，因为$F\equiv F_0\pmod{x^{\frac{n}{2}}}$，所以后面部分一定有长度为$\frac{n}{2}$的部分全部都是$0$。此时，我们可以将后半部分位移至前面，然后再进行操作。通常可以将多项式乘法的长度减少一半，时间复杂度会大大优化。
例如，多项式开方的公式为$F=F_0+\frac{1}{2}F_0^{-1}(G-F_0^2)$。应用以上方法可以知道，$(G-F_0^2)$有一半都是$0$。那么，如果当前层长度为$n$，从$\frac{n}{2}$转移过来，此时我们做多项式乘法只需要做长度为$n$的而不是$2n$的。同理，求逆也只需要处理到$n$。这样可以省下近一半常数。
## 8.多项式 $\exp$
[点我](https://www.luogu.org/problemnew/show/P4726)
由泰勒展开可得 $e^x=\sum_{i=0}^\infty \frac{x^i}{i!}$。因此，我们定义 $\sum_{i=0}^\infty \frac{F(x)^i}{i!}=\exp{F}$。介于 $e$ 在膜域下无定义，一般默认$F(x)$ 常数项为 $0$ 。
已知函数 $F$ ，求 $G\equiv \exp{F}$ 。保证 $F$ 常数项为 $0$ 。
可以知道 $\exp(\ln P)=P$ 。
所以 $\ln G-F=0$ 。
直接套用上方多项式牛顿迭代的公式，可得
$G=G_0(1-\ln G_0+F)$
倍增即可。时间复杂度 $O(n\log n)$ 。
该算法常数巨大，经验算大约在 $70\sim 80$ 左右。因此可以粗略认为是 $O(n\log^2 n)$ 的。
~~最近发现常数过大还没有[分治 $FFT$](https://www.luogu.org/recordnew/show/17449262) 跑得快。~~
#### 代码:
```cpp
static int ExX[MAXN],ExY[MAXN];
inline void Exp(int *F,int *G,int ln)
{
G[0]=1;
for(register int Lx=2;Lx<=ln;Lx<<=1)
{
rep(i,Lx>>1,Lx)G[i]=0;
Ln(G,ExX,Lx);
rep(i,0,Lx>>1)ExX[i]=ad(F[i+(Lx>>1)],mod-ExX[i+(Lx>>1)]);
rep(i,0,Lx>>1)ExY[i]=G[i];
rep(i,Lx>>1,Lx)ExX[i]=0;
Len=Lx;
calrev();
NTT(ExY,1),NTT(ExX,1);
rep(i,0,Len)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
NTT(ExX,-1);
rep(i,0,Len>>1)G[i+(Len>>1)]=ExX[i];
}
}
```
## 9.多项式幂次
例题见10.拉格朗日反演。
直接上式子吧...
当 $F[0]=1$ 时:
$$F^k=\exp(k\ln F)$$
否则，设最低项为 $a_tx^t$ ,则
$$F^k=a_t^kx^{tk}\exp(k\ln \frac{F}{a_tx^t})$$
直接按照式子计算即可。时间复杂度$O(n\log n)$ 。
#### 注意:
直接求 $\ln$ 再 $\exp$ 的方法仅适用于常数项为1的情况(因为多项式 $\exp$ 计算时默认常数项为 $1$ )。因此当常数项不为 $1$ 时需要将多项式的常数项强制变成 $1$ 。多判一下即可(不过大部分时间都不用管)。
#### 代码:
```cpp
inline void Pow(int*F,int*G,int ln,ll k)
{
int lst=ln+1;
Rep(i,0,ln)if(F[i]){lst=i;break;}
if(ln&&(__int128)lst*k>ln)
{memset(G,0,sizeof(int)*(ln+1));return;}
int md=ln-k*lst,bs=F[lst],iv=power(bs,mod-2);
Rep(i,lst,ln)G[i]=(ll)G[i]*iv%mod;
Ln(G+lst,G,md);
k%=mod;
Rep(i,0,md)G[i]=(ll)G[i]*k%mod;
Exp(G,G,md);
bs=power(bs,k);
Repe(i,md,0)G[i+lst*k]=(ll)G[i]*bs%mod;
memset(G,0,sizeof(int)*(lst*k));
}
```
### 9.5.多项式 $k$ 次方根
在膜意义下，存在一下转化:
$$\sqrt[k]F=F^{\frac{1}{k}}=F^{k^{-1}}$$
直接套用上面的幂次即可。复杂度为 $O(n\log n)$ 。
代码就咕了。
## 10.拉格朗日反演
[点我](https://www.lydsy.com/JudgeOnline/problem.php?id=3684)
已知 $G$ ，求函数 $F$ 使得 $F[G(x)]=x$ 。
经过[复杂的推理](http://zjt-blog.cc/articles/63?jump_page=2)，存在关系式:
$$[x^n] F(x)=\frac{1}{n}[x^{n-1}](\frac{x}{G(x)})^n$$
其中$[x^n]$表示多项式$x^n$项系数。
直接按照式子模拟即可。时间复杂度 $O(n\log n)$ 。
#### 代码:
```cpp
inline int Lagrange(int *a,int n)
{
Rep(i,0,n-1)a[i]=a[i+1];
Inv(a,n);
pow(a,n,n);
return (ll)a[n-1]*power(n,mod-2)%mod;
}
```
### 10.5 扩展拉格朗日反演
其实就是另一个式子。
$$[x^n]H[F(x)]=\frac{1}{n}[x^{n-1}]H'(x)(\frac{x}{G(x)})^n$$
## 11.分治FFT
[点我](https://www.luogu.org/problemnew/show/P4721)
已知 $F[0]=1$ ,求 $F$ 满足
$$F[i]=\sum_{j=1}^i F[i-j]*G[j]$$
其实哪怕是$F[i]=aF[i-b]+c+\sum_{j=1}^i dF[i-j]*G[j]$也是可以解决的。
考虑利用分治。
假设我们求出了$l\to mid$的答案，要求这些点对$mid+1\to r$的影响，那么对右半边点$x$的贡献为:
$w_x=\displaystyle\sum_{i=l}^{mid} f[i]g[x-i]$
这部分可以利用卷积来快速计算。计算完以后，答案直接加到答案数组就可以了。
需要注意的是，如果要求左边点对右边点的影响，首先整个区间以左对该区间的贡献应该先求出。所以分治过程为先分治左边，在求出中间，然后在递归右边。
时间复杂度 $O(n\log^2n)$ 。
有一个简单卡常方法。求一层的时候，我们会做长度为 $md-l+1+r-l+1\approx 1.5(r-l+1)$ 的卷积，但是这没有必要。一个显然的观察是我们只会使用在 $md-l+1\sim r-l+1$ 内的值，而前半截内容显然是没有用的。因此我们可以使用长度为 $r-l+1$ 的卷积。一般情况下可以快 $50\%$ 。
代码(求 $\exp$):
```cpp
void cdq_FFT(int l,int r)
{
if(l==r)
{
if(!l)F[l]=1;
else F[l]=(ll)F[l]*iv[l]%mod;
return;
}
int md=(l+r)>>1;cdq_FFT(l,md);
for(Len=2;Len<=r-l;Len<<=1);
calrev();
memset(X,0,sizeof(int)*Len);
memset(Y,0,sizeof(int)*Len);
memcpy(X,F+l,sizeof(int)*(md-l+1));
memcpy(Y,a,sizeof(int)*(r-l));
NTT(X,1),NTT(Y,1);
rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
NTT(X,-1);
Rep(i,md+1,r)F[i]=ad(F[i],X[i-l-1]);
cdq_FFT(md+1,r);
}
```
### 11.5.特殊的分治FFT
[点我](https://www.luogu.org/problemnew/show/P4566)
先假设你会前面的推式子部分。
你现在得到了一个长得很丑的式子，大概长这样:
$$f_n=(n-1)f_{n-1}+\sum_{i=1}^{n-2}(i-1)f_if_{n-i}$$
换句话说，你得到了一个类似于$F=F+F*F$ 的式子。这个式子利用上面的分治 $FFT$ 方法当然是求不出的了。
这时候，你就需要对你的分治 $FFT$ 进行魔改了。魔改结果大概长这样:
```cpp
void cdqFFT(int *f,int l,int r)
{
if(l==r)
{
f[l]=ad(f[l],(uint64)(l-1)*f[l-1]%mod);
return;
}
int mid=(l+r)>>1;
cdqFFT(f,l,mid);
for(Len=r-l+1;Len>(Len&-Len);Len+=(Len&-Len));
Rep(i,0,Len)A[i]=B[i]=0;
Rep(i,l,mid)A[i-l]=(uint64)f[i]*(i-1)%mod,B[i-l]=f[i];
calrev(Len);mul(A,B);
Rep(i,mid+1,r)if(i>=l*2)f[i]=ad(f[i],A[i-2*l]);
if(l^2)
{
Rep(i,0,Len)A[i]=B[i]=0;
Rep(i,2,min(l-1,r-l))A[i-2]=f[i];
Rep(i,l,mid)B[i-l]=f[i];
mul(A,B);
Rep(i,mid+1,r)if(i>=l+2)
f[i]=ad(f[i],(uint64)(i-2)*A[i-l-2]%mod);
}
cdqFFT(f,mid+1,r);
}
```
稍微解释一下。
首先，前面的 $(n-1)f_{n-1}$ 明显是可以在递归最底层在加上的。这个很好处理。
问题就在于后面的 $\sum_{i=1}^{n-2}(i-1)f_if_{n-i}$。这部分很难求。但是你可以分两种情况来进行讨论:
1. 所计算区间左端点是全局左端点。原来怎么做就怎么做。但是注意不要在卷积中用到还没有计算完毕的项。
2. 所计算区间不是全局左端点。可以知道会存在一些项之前并没有计算过。我们将这些项合并一起计算。具体来说就是把 $(i-1)f_if_{n-i}$ 和 $(n-i-1)f_{n-i}f_i$ 一起计算。那么这两项加在一起就变成了$(n-2)f_if_{n-i}$ 。直接加上即可。
其余部分和一般分治$FFT$相同。注意处理几个边界情况。
注意，我程序内部有一个 $if(l\otimes 2)$ ，这是判断左端点是否为全局左端点的。因为这道题的全局左端点是 $2$ 。
## 12.常系数齐次线性递推
[点我](https://www.luogu.org/problemnew/show/P4723)
主要就是推式子。
推导过程和代码见[线性递推](https://www.luogu.org/blog/user7035/solution-p4723)。
## 13.多项式求三角函数
经过审稿管理员提示，我发现原来这东西还是[有用](https://en.wikipedia.org/wiki/Alternating_permutation)的。。。
已知$F(x)$，求$p[F(x)]$。其中， $p(x)$ 为某一个三角函数 $(\sin,\cos,\tan,\sec,\csc,\cot)$ 。
首先，可以知道，如果求出了$\sin$ 和$\cos$ ，那么其它的三角函数也可以轻易推出。那么问题就是如何求出这两个三角函数。
根据欧拉定理可知:
$$e^{w_4x}=\cos(x)+w_4\sin(x)$$
其中， $w_4$ 为四次单位根，也就是常用的 $i$ 。
那么同理可得
$$e^{w_4F(x)}=\cos[F(x)]+w_4\sin[F(x)]$$
将$F(x)$ 的系数都乘以 $w_4$ 然后 $\exp$ 就可以了。时间复杂度 $O(n\log n)$ 。
注意以上的内容是 $FFT$ 的求法。如果要用 $NTT$ 的话也不是不行。
由
$$\cos(x)=\sum_{i=0}^\infty \frac{[i\bmod 2=0](-1)^\frac{i}{2}x^i}{i!}$$
$$\cosh(x)=\frac{e^x+e^{-x}}{2}=\sum_{i=0}^\infty \frac{[i\bmod 2=0]x^i}{i!}$$
可得
$$\cos(x)=\frac{e^{w_4x}+e^{-w_4x}}{2}$$
$$\cos(F)=\frac{1}{2}(e^{w_4F}-e^{-w_4F})$$
因为 $w_4$ 在常用 $NTT$ 模数下一般都有定义，且实系数函数的 $\cos$ 肯定也是实系数函数，因此这个公式计算出来肯定不含复系数。
至于 $\sin$ ，可以利用类似的方法来计算。即:
$$\sin(x)=\sum_{i=0}^\infty \frac{[i\bmod2=1](-1)^{\frac{i-1}{2}}x^i}{i!}$$
$$\sinh(x)=\frac{e^x-e^{-x}}{2}=\sum_{i=0}^\infty \frac{[i\bmod2=1]x^i}{i!}$$
推得:
$$\sin(x)=\frac{w_4^3}{2}(e^{w_4x}-e^{-w_4x})$$
同理得:
$$\sin F=\frac{w_4^3}{2}(e^{w_4F}-e^{-w_4F})$$
#### 代码:
```cpp
static int Sn[MAXN],Cs[MAXN];
inline void Cos(int*F,int*G,int ln)
{
Rep(i,0,ln)Cs[i]=(ll)F[i]*w4%mod;
Exp(Cs,Cs,ln),Inv(Cs,G,ln);
Rep(i,0,ln)G[i]=((ll)Cs[i]+G[i])*inv2%mod;
}
const int w43div2=(ll)w4*w4%mod*w4%mod*inv2%mod;
inline void Sin(int*F,int*G,int ln)
{
Rep(i,0,ln)Sn[i]=(ll)F[i]*w4%mod;
Exp(Sn,Sn,ln),Inv(Sn,G,ln);
Rep(i,0,ln)G[i]=((ll)Sn[i]+mod-G[i])*w43div2%mod;
}
```
## 14.多项式反三角函数
这个东西和上面的没有半毛钱关系。。。
根据积分表可以知道:
$$\arcsin x = \int \frac{1}{\sqrt{1-x^2}}$$
$$\arctan x = \int \frac{1}{1+x^2}$$
那么同理可得:
$$\arcsin F = \int \frac{F'}{\sqrt{1-F^2}}$$
$$\arctan F = \int \frac{F'}{1+F^2}$$
直接计算即可。
代码:
```cpp
inline void Arcsin(int*F,int*G,int ln)
{
memcpy(Sn,F,sizeof(int)*(ln+1));
for(Len=2;Len<=ln*2;Len<<=1);calrev();
NTT(Sn,1);rep(i,0,Len)Sn[i]=(ll)Sn[i]*Sn[i]%mod;
NTT(Sn,-1);
Sn[0]=1;rep(i,ln+1,Len)Sn[i]=0;
Rep(i,1,ln)Sn[i]=mod-Sn[i];
Sqrt(Sn,Sn,ln);
Inv(Sn,Sn,ln);
Deriv(F,G,ln);
mul(G,Sn,G,ln,ln);
Rep(i,ln+1,ln<<1)G[i]=0;
Inter(G,G,ln);
}
inline void Arctan(int*F,int*G,int ln)
{
memcpy(Sn,F,sizeof(int)*(ln+1));
for(Len=2;Len<=ln*2;Len<<=1);calrev();
NTT(Sn,1);rep(i,0,Len)Sn[i]=(ll)Sn[i]*Sn[i]%mod;
NTT(Sn,-1),Sn[0]=1;
rep(i,ln+1,Len)Sn[i]=0;Inv(Sn,Sn,ln);
Deriv(F,G,ln),mul(G,Sn,G,ln,ln);
Rep(i,ln+1,ln<<1)G[i]=0;
Inter(G,G,ln);
}
```
## 15.多项式多点求值
[点我](https://www.luogu.org/problemnew/show/P5050)
已知函数 $F$ 和数组 $\{a_i\}$ ，求 $\{F(a_i)\}$ 。
直接计算是 $O(nm)$ 的。此时，我们可以考虑降幂。
构造一个函数 $D_{l,r}=\prod_{i=l}^r (x-a_i)$ 。
可以发现，$\forall k\in [l,r],D_{l,r}(a_k)=0 $。
因此，$\forall k\in [l,r],F(a_k)=F(a_k)\bmod D_{l,r}(a_k) $。
于是，设 $F_{l,r}(x)$ 表示 $a_l$ 到 $a_r$ 在原函数上的点构成的多项式，则:
$$F_{l,mid}(x)=F_{l,r}(x)\bmod D_{l,mid}(x)$$
直接分治即可。时间复杂度为 $O(n\log^2 n)$ 。注意 $D$ 的计算可以直接分治乘法每次乘起来，不影响复杂度。
#### 代码:
```cpp
static int solv[18][2][MAXN];
inline void calc(int*a,int l,int r,int lev,int dir)
{
if(l==r)
{
solv[lev][dir][0]=mod-a[l];
solv[lev][dir][1]=1;return;
}
int md=(l+r)>>1;
calc(a,l,md,lev+1,0),calc(a,md+1,r,lev+1,1);
mul(solv[lev+1][0],solv[lev+1][1],solv[lev][dir],md-l+1,r-md);
}
static int pol[18][MAXN],Q[MAXN];
void getnum(int*a,int*ans,int l,int r,int lev,int n)
{
if(n>r-l)
{
calc(a,l,r,0,0);
Div(pol[lev-(lev>0)],solv[0][0],Q,pol[lev],n,r-l+1);
n=r-l;
}
else if(lev){Rep(i,0,n)pol[lev][i]=pol[lev-1][i];}
if(l==r)
{
ans[l]=pol[lev][0];
return;
}
int md=(l+r)>>1;
getnum(a,ans,l,md,lev+1,n);
getnum(a,ans,md+1,r,lev+1,n);
}
```
## 16.多项式快速插值
[点我](https://www.luogu.org/problemnew/show/P5158)
给定二维平面上点集$\{(x_i,y_i)\}$(保证$x_i$互不相同)，求函数$F(x)$，满足所有点均在该函数上。
核心内容是优化拉格朗日插值。
具体过程见[多项式快速插值](https://www.luogu.org/blog/user7035/solution-p5158)。
## Ex1.MTT
[点我](https://www.luogu.org/problemnew/show/P4245)还有[我](https://www.luogu.org/problemnew/show/P4239)
上面的所有函数基于模数为一个 $NTT$ 模数(即 $a*2^k+1$ ，其中$k$比较大)。如果不是的话，可以将多项式乘法换成以下的函数:
```cpp
inline void mul(int*F,int*G)
{
rep(i,0,Len)A[i]=comp(F[i]>>15,F[i]&32767),
C[i]=comp(G[i]>>15,G[i]&32767);
DBLFFT(A,B,1),DBLFFT(C,D,1);
rep(i,0,Len)
{
register comp t=A[i];
A[i]=B[i]*C[i]+A[i]*D[i];
B[i]=B[i]*D[i];
C[i]=C[i]*t;
}
DBLFFT(A,B,-1),FFT(C,-1);
rep(i,0,Len)
{
register int lz=(((ll)floor(C[i].re+0.5))%mod*73741817
+((ll)floor(A[i].re+0.5))%mod*32768+(ll)floor(B[i].re+0.5))%mod;
F[i]=i<=k?lz:0;
}
}
```
其中 $DBLFFT$ 是共轭优化过的 $FFT$ 。下面一章有代码。
简单解释其实就是将系数拆成 $2^{15}a+b$ ，然后利用等式 $(2^{15}A+B)(2^{15}C+D)=2^{30}AC+2^{15}(BC+AD)+BD$ 直接计算。原理为 $double$ 的精度可以承受住 $10^{14}$ ，因此能够支撑 $2^{15}*2^{15}*10^5=107374182400000$ 以内的数字运算。
具体还是看 $2016$ 毛爷爷写的国家队论文吧。
顺便一提，配合共轭优化后可以做到 $4$ 次 $DFT$ ，基本和一般的乘法没什么区别了。
## Ex2.BlueStein算法
[这个算例题吧](https://www.luogu.org/problemnew/show/P5293)
已知多项式 $F(x)$ ，求 $\forall i\in[0,k),A_i=F(w_k^i)$ 。其中 $w_k$ 为 $k$ 次单位根。
被这个东西送到了退役边缘。。。只怪我自己没有认真看论文。。。
因为 $w_k^t=w_k^{t\bmod k}$ ，因此可以认为 $F$ 的最高次幂为 $k-1$ 。
设 $F$ 的 $x^t$ 次方系数为 $f_t$ 。
那么，将式子写出，得:
$$F(w_k^t)=\sum_{i=0}^{k-1}f_iw_k^{ti}$$
$$=\sum_{i=0}^{k-1}f_iw_k^{-(k-t)i}$$
设这个东西为 $ans_{k-t}$ ，则 $F(w_k^t)=ans_{k-t}$ 。
运用套路，得:
$$ti=\binom{t+i}{2}-\binom{t}{2}-\binom{i}{2}$$
因此有:
$$ans_t=\sum_{i=0}^{k-1} f_iw_k^{\binom{i-t}{2}-\binom{-t}{2}-\binom{i}{2}}$$
$$=w_k^{-\binom{-t}{2}}\sum_{i=0}^{k-1} f_iw_k^{-\binom{i}{2}}*w_k^{\binom{-(t-i)}{2}}$$
因为我们每一项都要卷到 $k-1$ ，因此考虑将系数直接后移 $k-1$ 位，可以得到新的式子:
$$ans_t=w_k^{-\binom{-t}{2}}\sum_{i=0}^{k-1+t} f_iw_k^{-\binom{i}{2}}*w_k^{\binom{k-1-(k-1+t-i)}{2}}$$
我们设 $A_i=f_iw_k^{-\binom{i}{2}}$ ， $B_i=w_k^{\binom{k-1-(k-1+t-i)}{2}}$ ，则
$$ans_t=w_k^{-\binom{-t}{2}}[k-1+t](A*B)$$
时间复杂度 $O(n\log n)$ 。
这个过程被称为 $CZT$ 。 $ICZT$ 只需要将 $F_1$ 到 $F_{k-1}$ 的部分翻转，然后整体除掉 $k$ 即可。
算出这个东西后，我们就可以将他当做以 $w_k$ 为公比的插值。因此如果将多个式子相卷，得到的结果将会是这些式子在膜 $x^k$ 意义下的循环卷积。有的时候会有用。
代码:
```cpp
inline void Blue_Stein(int*F,int*G,int ln)
{
rep(i,0,ln)A[i]=(ll)F[i]*power(ome,ln-(ll)i*(i-1)/2%ln,mod)%mod;
Rep(i,0,2*ln-2)B[i]=power(ome,ln+(ll)(ln-1-i)*(ln-2-i)/2%ln,mod);
mul(A,B,A,ln-1,2*ln-2);
reverse(A+1,A+k);
rep(i,0,ln)G[i]=(ll)A[i+k-1]*power(ome,ln-(ll)(-i)*(-i-1)/2%ln,mod)%mod;
}
```
注意使用的时候大部分情况下模数都不会很友好，因此结合上面的 $MTT$ 或者直接上 $FFT$ 是一般情况。
## Ex3.常数优化
根据经验可以发现，多项式方面的常数差异非常大，有的时候别人跑过了但是你没跑过就凉了。因此，常数优化是非常重要的。限于篇幅就只推荐[这篇](http://negiizhao.blog.uoj.ac/blog/4671)和[这篇](http://zjt-blog.cc/articles/154?jump_page=2)了。后面那篇的详细内容和不详细证明可以参考毛爷爷的国家队论文。
此处放出共轭优化的代码。该代码可以同时对 $F$ 和 $G$ 进行 $DFT$ 或者 $IDFT$ 操作。
```cpp
static comp Z[MAXN];
inline void DBLFFT(comp*F,comp*G,int tl)
{
if(tl==1)
{
memcpy(Z,F,sizeof(comp)*Len);
FFT(Z,1);
memcpy(F,Z,sizeof(comp)*Len);
reverse(F+1,F+Len);
rep(i,0,Len)F[i].im=-F[i].im;
rep(i,0,Len)G[i]=(Z[i]-F[i])/2,F[i]=(Z[i]+F[i])/2
,swap(G[i].re,G[i].im),G[i].im=-G[i].im;
}
else
{
rep(i,0,Len)F[i]=comp(F[i].re-G[i].im,F[i].im+G[i].re);
FFT(F,-1);
rep(i,0,Len)G[i]=comp(F[i].im,0),F[i]=comp(F[i].re,0);
}
}
```
## Ex4.快速子集变换
[点我](https://www.luogu.org/problemnew/show/P4717)
和 $FFT$ 的形式差不多。代码和内容在[这里](https://www.luogu.org/blog/user7035/solution-p4717)。
## Ex4.5快速子集运算
限于篇幅在[另一篇文章](https://www.luogu.org/blog/user7035/zi-ji-juan-ji-ji-ji-gao-ji-yun-suan)作说明。
## 总集:[帕秋莉的超级多项式](http://cogs.pro:8080/cogs/problem/problem.php?pid=2189)
已知 $n-1$ 次多项式 $F$ 和 $k$ ，求
$$\left[\left(1+\ln(1+\frac{1}{\exp(\int\frac{1}{\sqrt{F}})})\right)^k\right]'$$
这道题汇集了几乎所有的多项式基础操作，有兴趣的可以尝试做一下，练练手。
#### 代码:
```cpp
#include<bits/stdc++.h>
#define Rep(i,a,b) for(register int i=(a);i<=(b);++i)
#define Repe(i,a,b) for(register int i=(a);i>=(b);--i)
#define rep(i,a,b) for(register int i=(a);i<(b);++i)
#define pb push_back
#define mx(a,b) (a>b?a:b)
#define mn(a,b) (a<b?a:b)
typedef unsigned long long uint64;
typedef unsigned int uint32;
typedef long long ll;
using namespace std;
namespace IO
{
const uint32 Buffsize=1<<15,Output=1<<24;
static char Ch[Buffsize],*S=Ch,*T=Ch;
inline char getc()
{
return((S==T)&&(T=(S=Ch)+fread(Ch,1,Buffsize,stdin),S==T)?0:*S++);
}
static char Out[Output],*nowps=Out;
inline void flush(){fwrite(Out,1,nowps-Out,stdout);nowps=Out;}
template<typename T>inline void read(T&x)
{
x=0;static char ch;T f=1;
for(ch=getc();!isdigit(ch);ch=getc())if(ch=='-')f=-1;
for(;isdigit(ch);ch=getc())x=x*10+(ch^48);
x*=f;
}
template<typename T>inline void write(T x,char ch='\n')
{
if(!x)*nowps++='0';
if(x<0)*nowps++='-',x=-x;
static uint32 sta[111],tp;
for(tp=0;x;x/=10)sta[++tp]=x%10;
for(;tp;*nowps++=sta[tp--]^48);
*nowps++=ch;
}
}
using namespace IO;
void file(void)
{
FILE*DSD=freopen("polynomial.in","r",stdin);
FILE*CSC=freopen("polynomial.out","w",stdout);
}
const int MAXN=1<<21;
namespace poly
{
const int mod=998244353,gen=3;
static int g[23],iv[MAXN];
inline int power(int u,int v)
{
register int sm=1;
for(;v;v>>=1,u=(uint64)u*u%mod)if(v&1)
sm=(uint64)sm*u%mod;
return sm;
}
inline void predone()
{
Rep(i,1,21)g[i]=power(gen,(mod-1)>>i);
iv[1]=1;
Rep(i,2,1e5)iv[i]=mod-(uint64)mod/i*iv[mod%i]%mod;
}
static int Len,rev[MAXN];
inline void calrev()
{
int II=log(Len)/log(2)-1;
Rep(i,1,Len-1)rev[i]=rev[i>>1]>>1|(i&1)<<II;
}
inline int ad(int u,int v){return(u+=v)>=mod?u-mod:u;}
inline void NTT(int*F,int typ)
{
Rep(i,1,Len-1)if(i<rev[i])swap(*(F+i),*(F+*(rev+i)));
for(register int i=2,ii=1,t=1;i<=Len;i<<=1,ii<<=1,++t)
{
register int wn=*(g+t);
for(register int j=0;j<Len;j+=i)
{
register int w=1;
rep(k,0,ii)
{
register int tt=(uint64)*(F+j+k+ii)*w%mod;
*(F+j+k+ii)=ad(*(F+j+k),mod-tt);
*(F+j+k)=ad(*(F+j+k),tt);
w=(uint64)w*wn%mod;
}
}
}
if(typ==-1)
{
reverse(F+1,F+Len);
register uint64 invn=power(Len,mod-2);
rep(i,0,Len)*(F+i)=invn**(F+i)%mod;
}
}
static int X[MAXN],Y[MAXN],Iv[MAXN];
inline void mul(int *a,int *b,int *c,int lenl,int lenr)
{
if((ll)lenl*lenr<=300)
{
Rep(i,0,lenl+lenr)X[i]=0;
Rep(i,0,lenl)Rep(j,0,lenr)
X[i+j]=(X[i+j]+(ll)a[i]*b[j])%mod;
Rep(i,0,lenl+lenr)c[i]=X[i];
return;
}
for(Len=2;Len<=lenl+lenr;Len<<=1);
calrev();
Rep(i,0,lenl)X[i]=a[i];
Rep(i,0,lenr)Y[i]=b[i];
rep(i,lenl+1,Len)X[i]=0;
rep(i,lenr+1,Len)Y[i]=0;
NTT(X,1),NTT(Y,1);
rep(i,0,Len)X[i]=(ll)X[i]*Y[i]%mod;
NTT(X,-1);
Rep(i,0,lenl+lenr)c[i]=X[i];
rep(i,lenl+lenr+1,Len)c[i]=0;
}
inline void Inv(int*F,int*G,int ln)
{
Iv[0]=power(F[0],mod-2);
for(register int Ln=2;Ln>>1<=ln;Ln<<=1)
{
rep(i,ln+1,Ln)F[i]=0;
rep(i,0,Ln)X[i]=F[i],Y[i]=0;
rep(i,0,(Ln>>1))Y[i]=Iv[i];
Len=Ln,calrev();
NTT(X,1),NTT(Y,1);
rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
NTT(X,-1);
rep(i,0,(Ln>>1))X[i]=0;
NTT(X,1);
rep(i,0,Ln)X[i]=(uint64)X[i]*Y[i]%mod;
NTT(X,-1);
rep(i,(Ln>>1),Ln)Iv[i]=mod-X[i];
}
Rep(i,0,ln)G[i]=Iv[i];
}
static int ExX[MAXN],ExY[MAXN],Op[MAXN];
inline void Div(int*F,int*G,int*Q,int*R,int lenf,int leng)
{
Rep(i,0,lenf)ExX[i]=F[lenf-i];
Rep(i,0,leng)ExY[i]=G[leng-i];
Rep(i,leng+1,lenf-leng)ExY[i]=0;
Inv(ExY,ExY,lenf-leng);
Rep(i,lenf-leng+1,lenf)ExX[i]=0;
Rep(i,lenf-leng+1,leng)ExY[i]=0;
mul(ExX,ExY,ExY,lenf-leng,lenf-leng);
Rep(i,lenf-leng+1,(lenf-leng)<<1)ExY[i]=0;
Rep(i,0,(lenf-leng)>>1)swap(ExY[i],ExY[lenf-leng-i]);
mul(ExY,G,ExX,lenf-leng,leng);
assert(F[leng]==ExX[leng]);
Rep(i,0,leng-1)R[i]=ad(F[i],mod-ExX[i]);
Rep(i,0,lenf-leng)Q[i]=ExY[i];
}
const int inv2=(mod+1)>>1;
inline void Sqrt(int*F,int*G,int ln)
{
Op[0]=sqrt(F[0]);
for(register int Ln=2;Ln>>1<=ln;Ln<<=1)
{
mul(Op,Op,ExX,Ln>>1,Ln>>1);
rep(i,0,Ln>>1)ExX[i]=ad(mod-ExX[i+(Ln>>1)],F[i+(Ln>>1)]);
rep(i,Ln>>1,Ln)Op[i]=ExX[i]=0;
Inv(Op,ExY,Ln>>1);
Len=Ln,calrev();
NTT(ExX,1),NTT(ExY,1);
rep(i,0,Ln)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
NTT(ExX,-1);
rep(i,0,Ln>>1)Op[i+(Ln>>1)]=(ll)ExX[i]*inv2%mod;
}
Rep(i,0,ln)G[i]=Op[i];
}
inline void Deriv(int*F,int*G,int ln)
{Rep(i,1,ln)G[i-1]=(uint64)F[i]*i%mod;G[ln]=0;}
inline void Inter(int*F,int*G,int ln)
{Repe(i,ln,1)G[i]=(uint64)F[i-1]*iv[i]%mod;G[0]=0;}
static int LnX[MAXN];
inline void Ln(int*F,int*G,int ln)
{
Deriv(F,LnX,ln),Inv(F,G,ln);
for(Len=2;Len<=ln<<1;Len<<=1);
rep(i,ln+1,Len)LnX[i]=G[i]=0;
calrev();
NTT(LnX,1),NTT(G,1);
rep(i,0,Len)G[i]=(uint64)G[i]*LnX[i]%mod;
NTT(G,-1);
Inter(G,G,ln);
}
inline void Exp(int *F,int *G,int ln)
{
Op[0]=1;
for(register int Lx=2;Lx>>1<=ln;Lx<<=1)
{
rep(i,Lx>>1,Lx)Op[i]=0;
Ln(Op,ExX,Lx);
rep(i,0,Lx>>1)ExX[i]=ad(F[i+(Lx>>1)],mod-ExX[i+(Lx>>1)]);
rep(i,0,Lx>>1)ExY[i]=Op[i];
rep(i,Lx>>1,Lx)ExX[i]=ExY[i]=0;
Len=Lx;
calrev();
NTT(ExY,1),NTT(ExX,1);
rep(i,0,Len)ExX[i]=(ll)ExX[i]*ExY[i]%mod;
NTT(ExX,-1);
rep(i,0,Lx>>1)Op[i+(Len>>1)]=ExX[i];
}
Rep(i,0,ln)G[i]=Op[i];
}
}
using poly::power;
using poly::Len;
using poly::calrev;
using poly::NTT;
using poly::mod;
using poly::predone;
using poly::Inv;
using poly::Inter;
using poly::Deriv;
using poly::Ln;
using poly::Exp;
using poly::Sqrt;
using poly::Div;
using poly::ad;
static int n,k,F[MAXN];
int main()
{
file();
predone();
read(n),read(k);
k%=mod;
Rep(i,0,n-1)read(F[i]);
Sqrt(F,F,n-1);
Inv(F,F,n-1);
Inter(F,F,n-1);
Exp(F,F,n-1);
Inv(F,F,n-1);
++F[0];
Ln(F,F,n-1);
++F[0];
Ln(F,F,n-1);
Rep(i,0,n-1)F[i]=(ll)F[i]*k%mod;
Exp(F,F,n-1);
Deriv(F,F,n-1);
Rep(i,0,n-1)write(F[i],' ');
flush();
return 0;
}
```
推荐这种东西写成一个 $namespace$ ，方便调试~~和复板子~~。