多项式从入门到入坟
参考
基本概念
多项式一般有三种表达方式：
-
系数表达式：即用一个系数数组 \(A\) 表示这个多项式。
即\(F(x)=\sum_{i=0}^n A_ix^i\) -
点值表达式：因为 \(n+1\) 个点确定一个多项式，所以可以用一个 \(n+1\) 个点表示一个多项式，显然点值表达式不唯一。
-
下降幂表达式，可以通过斯特林数与系数表达式相互转化。
卷积
当两个多项式（设系数数组为 \(F,G\) 分别为 \(n\) 次和 \(m\) 次），当二者相乘的时候，得到的多项式的第 \(i\) 项系数为：\(\sum_{j+k=i}F_j \times G_k\)
该形式就是卷积。
正常两个多项式相乘需要 \(O(n^2)\) 的复杂度，而利用多项式算法可以优化成 \(O(n\log_n)\) 。
FFT
接下来看如何优化卷积。
- 复数
基础概念请自行翻阅必修 2。
一个复数可以表示成 \(a+bi\) 的形式。
- 加法：\((a+bi)+(c+di)=(a+c)+(b+d)i\)
- 减法：\((a+bi)-(c+di)=(a-c)+(b-d)i\)
- 乘法：\((a+bi)\times(c+di)=(ac-bd)+(ad \times bc)i\)
两个复数相乘，模长相乘，幅角相加
struct cmplx{
double x,y;
cmplx(){}
cmplx(double a,double b){x=a,y=b;}
cmplx operator + (const cmplx &A)const{return cmplx(x+A.x,y+A.y);}
cmplx operator - (const cmplx &A)const{return cmplx(x-A.x,y-A.y);}
cmplx operator * (const cmplx &A)const{return cmplx(x*A.x-y*A.y,x*A.y+y*A.x);}
};
- 单位根
显然，\(x^n=1\) 在复数域内有且只有 \(n\) 个解，我们称这些解为 \(n\) 次单位根，其幅角为 \(\frac{a}{n}\times 2\pi(0 \le a <n,a \in Z)\)
记第 \(k\) 个单位根为 \(\omega_n^k\)。
单位根有以下性质(可以用几何意义来理解)：
\(\omega_n^i \times \omega_n^j=\omega_n^{i+j}\)
\(\omega_{pn}^{pk}=\omega_n^k\)
\(\omega_n^{k+n/2}=-\omega_n^k\)
\(\omega_n^k=\omega_n^{k\bmod n}\)
而因为
\(\omega_n^1= cos\frac{2\pi}{n}+sin\frac{2\pi}{n}i\)
\(\omega_n^k=(\omega_n^1)^k\)
就可以得到所有单位根。
- DFT
用于快速把系数表达式转成点值表达式。
以下的 \(n\) 都满足 \(2^k=n,k \in N\)。
对于一个多项式 \(A(x)=\sum_{i=0}^{n-1}x^iF(i)\)，
\(A(x)=F(0)+xF(1)+x^2F(2)+...+x^{n-1}F(n-1)\)。
我们把它奇偶分开。
\(A1(x)=F(0)+xF(2)+...+x^{\frac{n}{2}-1}F(n-2)\)
\(A2(x)=F(1)+xF(3)+...+x^{\frac{n}{2}-1}F(n-1)\)
显然 \(A(x)=A1(x^2)+xA2(x^2)\)
设 \(k<\frac{n}{2}\)
把 \(\omega_n^k\) 带入
\(A(\omega_n^k)=A1((\omega_n^k)^2)+\omega_n^kA2((\omega_n^k)^2)\)
\(=A1(\omega_n^{2k})+\omega_n^kA2(\omega_n^{2k})\)
\(=A1(\omega_{n/2}^k)+\omega_n^kA2(\omega_{n/2}^k)\)
把 \(\omega_n^{k+n/2}\) 带入
\(A(\omega_n^{k+n/2})=A1((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}A2((\omega_n^{k+n/2})^2)\)
\(=A1(\omega_n^{2k+n})+\omega_n^{k+n/2}A2(\omega_n^{2k+n})\)
\(\because \omega_n^k=\omega_n^{k\ mod\ n}\)
\(\therefore A(\omega_n^{k+n/2})=A1(\omega_n^{2k})+\omega_n^{k+n/2}A2(\omega_n^{2k})\)
\(=A1(\omega_{n/2}^k)+\omega_n^{k+n/2}A2(\omega_{n/2}^k)\)
\(=A1(\omega_{n/2}^k)-\omega_n^kA2(\omega_{n/2}^k)\)
综上所述：
\(A(\omega_n^k)=A1(\omega_{n/2}^k)+\omega_n^kA2(\omega_{n/2}^k)\)
\(A(\omega_n^{k+n/2})=A1(\omega_{n/2}^k)-\omega_n^kA2(\omega_{n/2}^k)\)
就可以用分治在 \(O(n\log_n)\) 的时间内求出 \(n\) 个点值，就得到 \(A\) 的点值表达式。
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const int N=1e6+10;
const double pi=acos(-1);
struct cmplx
{
double x,y;
cmplx (double a=0,double b=0){ x=a,y=b; }
cmplx operator + (cmplx const &a) const{ return cmplx(x+a.x,y+a.y); }
cmplx operator - (cmplx const &a) const{ return cmplx(x-a.x,y-a.y); }
cmplx operator * (cmplx const &a) const{ return cmplx(x*a.x-y*a.y,x*a.y+y*a.x); }
}v[N],sav[N];
inline void DFT(cmplx *f,int len)
{
if(len==1) return;
cmplx *f1=f,*f2=f+len/2;//用指针，改f1/f2，f也会跟着改；改f，f1/f2也会改（它们用同样的地址）
for(int i=0;i<len;i++)//所以要开临时的数组
sav[i]=f[i];
for(int i=0;i<len/2;i++)
f1[i]=sav[i<<1],f2[i]=sav[i<<1|1];//奇偶分开
DFT(f1,len/2),DFT(f2,len/2);//分治
cmplx bas=cmplx(1,0),buf=cmplx(cos(2*pi/len),sin(2*pi/len));//算单位根
for(int i=0;i<len/2;i++,bas=bas*buf)
sav[i]=f1[i]+bas*f2[i],sav[i+len/2]=f1[i]-bas*f2[i];//套那两个式子
for(int i=0;i<len;i++)
f[i]=sav[i];
}
int main()
{
int n,m=1;
scanf("%d",&n);
for(int i=0;i<n;i++)
scanf("%lf",&v[i].x);//输入系数
while(m<n) m<<=1;//补全
DFT(v,m);
return 0;
}
注意：\(n=2^k\)，\(n\) 为项数！
- IDFT
其实是DFT的逆运算。
设 \(F(x)\) DFT 之后得到 \(G(x)\)。
单位根反演乱搞可得：
长得很像 DFT，所以代码也差不多（只差一个符号）。
- FFT
因为DFT和IDFT长得太像，所以把他们合并成FFT。
而且三角函数可以预处理。
通过找规律，初始序列最后变换到的位置是可以算出来的,这就叫蝴蝶变换：
for(int i=0;i<len;i++) cir[i]=(cir[i>>1]>>1)|(i&1?len>>1:0);
所以我们就得到了FFT,下面是迭代版的。
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const int N=4e6+10;
int cir[N];
double tsin[N],tcos[N];
const double pi=acos(-1);
struct cmplx
{
double x,y;
cmplx (double a=0,double b=0){ x=a,y=b; }
cmplx operator + (cmplx const &a) const{ return cmplx(x+a.x,y+a.y); }
cmplx operator - (cmplx const &a) const{ return cmplx(x-a.x,y-a.y); }
cmplx operator * (cmplx const &a) const{ return cmplx(x*a.x-y*a.y,x*a.y+y*a.x); }
}F[N],G[N];
inline void FFT(cmplx *f,int len,double tag)
{
for(register int i=0;i<len;i++) if(i<cir[i]) swap(f[i],f[cir[i]]);
for(register int l=2;l<=len;l<<=1)
{
int size=l>>1;
cmplx buf=cmplx(tcos[l],tag*tsin[l]);
for(register int i=0;i<len;i+=l)
{
cmplx bas=cmplx(1,0);
for(register int j=i;j<i+size;j++,bas=bas*buf)
{
cmplx tmp=bas*f[j+size];
f[j+size]=f[j]-tmp;
f[j]=f[j]+tmp;
}
}
}
}
int main()
{
int n,m,len=1;
scanf("%d%d",&n,&m);
for(register int i=0;i<=n;i++)
scanf("%lf",&F[i].x);
for(register int i=0;i<=m;i++)
scanf("%lf",&G[i].x);
while(len<=n+m) len<<=1;
for(register int i=len;i>=1;i>>=1) tsin[i]=sin(2*pi/i),tcos[i]=cos(2*pi/i);
for(register int i=0;i<len;i++) cir[i]=(cir[i>>1]>>1)|((i&1)?len>>1:0);
FFT(F,len,1),FFT(G,len,1);
for(register int i=0;i<len;i++) F[i]=F[i]*G[i];
FFT(F,len,-1);
for(register int i=0;i<=n+m;i++) printf("%d ",(int)(F[i].x/len+0.49));
return 0;
}
- NTT
在 OI 中我们讨厌实数的运算，于是有了 NTT。
我们首先要找到一个东西替代单位根，于是找到了原根。
其实我们就是要找到一个东西像单位根一样转 \(n\) 次就回来的东西，这不就是原根吗，转 \(p-1\) 回来，那我们设置转一次是 \(\frac{p-1}{n}\)，那不就刚好转 \(n\) 次就回来了吗。
这需要 \(n\mid p-1\)，而 \(n=2^k\)，所以 \(p-1\) 中需要很多的 \(2\)。
一般我们都是用 \(998244353\) 为模数：
- \(998244353\) 大小合适，且是质数
- \(998244353-1=2^{23}\times 119\)
#include<iostream>
#include<stdio.h>
#define p 998244353
#define g 3
#define invg 332748118
#define N 4000005
#define int long long
using namespace std;
inline int read(){
int x=0,f=1;
char ch=getchar();
while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
return x*f;
}
int n,m,F[N],G[N],cir[N];
int power(int x,int b){
int res=1;
while(b){
if(b&1) res=res*x%p;
b>>=1;
x=x*x%p;
}
return res;
}
void NTT(int *f,int lenth,int tag){
for(int i=0;i<lenth;++i) if(i<cir[i]) swap(f[i],f[cir[i]]);
for(int len=2;len<=lenth;len<<=1){
int sze=len>>1,buf=power((tag?g:invg),(p-1)/len);
for(int i=0;i<lenth;i+=len){
int bas=1;
for(int j=i;j<i+sze;j++,bas=bas*buf%p){
int tmp=bas*f[j+sze]%p;
f[j+sze]=(f[j]-tmp+p)%p;
f[j]=(f[j]+tmp)%p;
}
}
}
}
signed main(){
n=read(),m=read();
for(int i=0;i<=n;++i) F[i]=read();
for(int i=0;i<=m;++i) G[i]=read();
int len=1;
while(len<=n+m) len*=2;
for(int i=0;i<len;++i) cir[i]=(cir[i>>1]>>1)|((i&1)?len>>1:0);
NTT(F,len,1),NTT(G,len,1);
for(int i=0;i<len;++i) F[i]=F[i]*G[i]%p;
NTT(F,len,0);
for(int i=0;i<n+m+1;++i) printf("%lld ",F[i]*power(len,p-2)%p);
return 0;
}
封装后的版本
#define mo 998244353
#define g 3
#define invg 332748118
int cir[N],rt;
int power(int x,int b){
int res=1;
while(b){
if(b&1) res=res*x%mo;
b>>=1;
x=x*x%mo;
}
return res;
}
void clear(int *A,int x,int y){
for(int i=x;i<y;++i) A[i]=0;
}
void copy(int *A,int *B,int x,int y){
for(int i=x;i<y;++i) A[i]=B[i];
}
void mul(int *A,int *B,int x,int y){
for(int i=x;i<y;++i) A[i]=A[i]*B[i]%mo;
}
void NTT(int *f,int lenth,int tag){
if(lenth!=rt)for(int i=0;i<lenth;++i)cir[i]=(cir[i>>1]>>1)|((i&1)?lenth>>1:0);
rt=lenth;
for(int i=0;i<lenth;++i) if(i<cir[i]) swap(f[i],f[cir[i]]);
for(int len=2;len<=lenth;len<<=1){
int sze=len>>1,buf=power((tag?g:invg),(mo-1)/len);
for(int i=0;i<lenth;i+=len){
int bas=1;
for(int j=i;j<i+sze;j++,bas=bas*buf%mo){
int tmp=bas*f[j+sze]%mo;
f[j+sze]=(f[j]-tmp+mo)%mo;
f[j]=(f[j]+tmp)%mo;
}
}
}
if(tag==1) return;
int invn=power(lenth,mo-2);
for(int i=0;i<lenth;++i) f[i]=f[i]*invn%mo;
}
int tmul[N];
void Polymul(int *A,int *B,int len1,int len2){
#define sav tmul
int len=1;
while(len<(len1<<1)) len<<=1;
copy(sav,B,0,len);
NTT(A,len,1),NTT(sav,len,1);
for(int i=0;i<len;++i) A[i]=A[i]*sav[i]%mo;
NTT(A,len,0);
clear(A,len2,len),clear(sav,0,len);
#undef sav
}
多项式求逆
设 \(G(x)\) 为答案，当 \(F(x)\) 只有1项时，答案为 \(F(x)\) 的逆元，其他情况考虑倍增。
设 \(mod\ x^{\lceil \frac{n}{2}\rceil}\) 时的答案为 \(H(x)\)
两边同乘 \(F(x)\)
int inv1[N],inv2[N],inv3[N];
void Polyinv(int *A,int len){
#define sav1 inv1
#define sav2 inv2
#define sav3 inv3
int lim=1;while(lim<len) lim<<=1;
sav1[0]=power(A[0],mo-2);
for(int l=2;l<=lim;l<<=1){
int r=l<<1;copy(sav3,A,0,l);
for(int i=0;i<(l>>1);++i) sav2[i]=(sav1[i]<<1)%mo;
NTT(sav1,r,1),mul(sav1,sav1,0,r);
NTT(sav3,r,1),mul(sav1,sav3,0,r);
NTT(sav1,r,0),clear(sav1,l,r);
for(int i=0;i<l;++i) sav1[i]=(sav2[i]-sav1[i]+mo)%mo;
}
copy(A,sav1,0,len);
clear(sav1,0,lim<<1),clear(sav2,0,lim<<1),clear(sav3,0,lim<<1);
#undef sav1
#undef sav2
#undef sav3
}
分治FFT
这题看似是 \(F=F*G\)。
但仔细想想却不是，因为如果按照这样 \(G[0]\) 应该是 \(0\) ，但题目又说 \(F[0]=1\)，所以并不能直接卷。
所以我们把 \(F 和 G\) 卷一下。
当 \(k=0\) ，那一项为\(0（G[0]=0）\)
当 \(k>0\) ，为 \(\sum_{i=0}^kF[k-i]G[i]=\sum_{i=1}^kF[k-i]G[i]=F[k]\) 。
所以 \(F*G\) 和 \(F\) 只有在 \(k=0\) 时差了个1，所以 \(F=F*G+1\) 。
移个项得： \(F=\dfrac{1}{1-G}\)，写个求逆就行了。
但是分治 FFT 的板子怎么能写求逆呢？而且求逆不是能完全代替分治 FFT 的。
分治 FFT 的本质是解方程，是一种半在线卷积。
考虑 \(F=F*G+1\)，模仿 CDQ 分治的过程，先算左边，再算左边对右边的贡献，最后再算右边。
具体地说，当我们在处理 \([l,r]\)，先 solve(l,mid)，然后把 \(F[l,mid]\) 和 \(G[0,r-l]\) 卷起来贡献到 \(F[mid+1,r]\)，再 solve(mid+1,r)。
但是有些题会出现自己卷自己的情况，例如 #50. 【UR #3】链式反应。
这题如果不管一些边界条件和细节的话，其实就是 \(F=\frac{1}{2}AF^2\)。
首先这不能求逆解决了，考虑分治 FFT，分两种情况考虑：
- \(l\not=0\)，这时候 \([0,r-l]\) 肯定已经算出来了，直接贡献即可。
- \(l=0\)，我们可以先只把 \([l,mid]\) 卷 \([l,mid]\) 的贡献算上，剩下的留给情况 1。所以情况 1 时不仅要考虑 \([l,mid]\) 对 \([mid+1,r]\) 的贡献，还要考虑 \([0,r-l]\) 对 \([mid+1,r]\) 的贡献，但因为 \([l,mid]\) 卷 \([0,r-l]\) 和 \([-,r-l]\) 卷 \([l,mid]\) 是一样的，所以我们需要做的只是在情况 1 时的贡献乘一个 2 而已。
void solve(int l,int r){
if(l==r) return;
int mid=l+r>>1;
solve(l,mid);
if(l!=0){
for(int i=l;i<=mid;++i) tmp1[i-l]=F[i];
for(int i=0;i<=r-l;++i) tmp2[i]=F[i];
polyMul(tmp1,tmp2,mid-l+1,r-l+1,r-l+1);
polyMul(tmp1,P,r-l+1,r-l+1,r-l+1);
for(int i=mid+1;i<=r;++i) (F[i]+=tmp1[i-l-1]*fac[i-1]%mo*ifac[i]%mo)%=mo;
clear(tmp1,0,r-l+1),clear(tmp2,0,r-l+1);
}
else{
for(int i=l;i<=mid;++i) tmp1[i-l]=F[i];
for(int i=l;i<=mid;++i) tmp2[i-l]=F[i];
polyMul(tmp1,tmp2,mid-l+1,mid-l+1,r-l+1);
polyMul(tmp1,P,r-l+1,r-l+1,r-l+1);
for(int i=mid+1;i<=r;++i) (F[i]+=tmp1[i-l-1]*fac[i-1]%mo*ifac[i]%mo*inv2%mo)%=mo;
clear(tmp1,0,r-l+1),clear(tmp2,0,r-l+1);
}
solve(mid+1,r);
}
其实形如 \(F=A*F\)，其中 \(A\) 是已知的，我们可以认为这是一次的类型。形如 \(F=A*F^2\)，这种是二次的类型，以此类推，其实二次就是自己卷自己的类型，所以上文讲了一次和二次的解法。
但方程不仅有几次，还有几元，比如 #428. 【集训队作业2018】普通的计数题。这题本质是一个二元二次：\(F=A*G,G=F*G\)。多元和一元的解法没多大区别，但是要注意，这里 \(l\not=0\) 时的贡献不能直接乘 2 了，因为原来我们直接乘 2 是因为 \(F[l,mid]*F[0,r-l]=F[0,r-l]*F[l,mid]\)，而这里 \(F[l,mid]*G[0,r-l]=G[l,mid]*F[0,r-l]\) 并不成立！所以我们只能重新算一遍，而不能直接乘以 2！
至于三次的解法其实很简单，例如 \(F=F^3\)，设 \(F_2=F^2,F=F_2*F\) 当二元二次解就好了。但是例如 \(F=F^8\)，其实不需要分别记 \(F_2,F_3\ldots F_8\)，其实只需要记 \(F_2,F_4,F_8\) 就好了，利用二进制拆分的原理可以优化一下常数。
多项式除法
给定 \(n\) 次多项式 \(F(x)\) 和 \(m\) 次的 \(G(x)\),求 \(n-m\) 次的 \(Q(x)\) 和 小于 \(m\) 次的 \(R(x)\)。
满足 \(F(x)\equiv G(x)Q(x)+R(x)\ (mod\ x^n)\)
思路是想办法让 \(R(x)\) 被模掉，但只有高次的会被模掉，而 \(R(x)\) 是最低次的。但是我们可以把所有多项式都反转，这样 \(R(x)\) 就是最高次的了！
设 \(F_0(x)\) 是 \(F(x)\) 反转之后的多项式。
显然 \(F_0(x)=\sum_{i=0}^{n}F[n-i]x^i\) 。
但我们需要更好的表达方式，发现 \(F_0(x)=x^nF(\dfrac{1}{x})\)，展开一下即可证明。
把上面的式子用 \(\dfrac{1}{x}\) 代替，并两边乘 \(x^n\) 。
左边就是 \(F_0(x)\)，而 \(G(x)\)是m次，\(Q(x)\)是n-m次，加起来刚好 n 次，\(R(x)\) 可以当成 m-1 次的，所以：
这下 \(x^{n-m+1}R_0(x)\) 是最高次了，所以：
求出 \(Q(x)\) 后代回去得出 \(R(x)\) 。
int div1[N],div2[N],rr[N];
void rev(int *A,int len){
#define sav rr
copy(sav,A,0,len);
for(int i=0;i<len;++i) A[i]=sav[len-i-1];
#undef sav
}
void Polydiv(int *A,int *B,int len1,int len2){
#define sav1 div1
#define sav2 div2
int len0=len1-len2+1;
rev(A,len1),copy(sav1,A,0,len0),rev(A,len1);
rev(B,len2),copy(sav2,B,0,len0),rev(B,len2);
Polyinv(sav2,len0),Polymul(sav2,sav1,len0,len0),rev(sav2,len0);
Polymul(B,sav2,len1,len1);
for(int i=0;i<len2-1;++i) B[i]=(A[i]-B[i]+mo)%mo;
copy(A,sav2,0,len0),clear(A,len0,len1),B[len2-1]=0;
clear(sav1,0,len1),clear(sav2,0,len1);
#undef sav1
#undef sav2
}
多项式求导与积分
不多说，大家都会。
int invcnt=1,inv[N];
inline void get_inv(int len){
inv[1]=1; if(invcnt>=len) return;
for(int i=invcnt+1;i<=len;++i) inv[i]=(mo-mo/i)*inv[mo%i]%mo;
}
inline void Polyder(int *A,int len){
for(int i=1;i<=len;i++) A[i-1]=A[i]*i%mo; A[len]=0;
}
inline void Polyint(int *A,int len){
get_inv(len);
for(int i=len;i>=1;i--) A[i]=A[i-1]*inv[i]%mo; A[0]=0;
}
多项式对数函数
给定 \(F(x)\) ，求 \(G(x)\) 使得 \(G(x) \equiv \ln F(x)\ \pmod {x^n}\)
因为 \(\ln\) 求导后很好看，所以两边求导。
int ln1[N];
inline void Polyln(int *A,int len){
#define sav ln1
copy(sav,A,0,len),Polyder(A,len),Polyinv(sav,len);
Polymul(A,sav,len,len),Polyint(A,len-1),clear(sav,0,len);
#undef sav
}
多项式牛顿迭代
求 \(F(x)\) 满足 \(G(F(x)) \equiv 0\ (mod\ x^n)\)
考虑递归求解，假设我们已经求出 \(\mod x^{\lceil\frac{n}{2}\rceil}\) 的解 \(F_0(x)\) 。
在 \(F_0(x)\) 处套泰勒展开：
已知 \(F(x)\) 和 \(F_0(x)\) 的后 \(\lceil\dfrac{n}{2}\rceil\) 项是一样的，所以 \(F(x)-F_0(x)\) 的大于等于 2 的次方不为零的项已经在 n 次以上了，在\(\mod x^n\) 情况下就是 0 ，所以：
移个项得：
多项式指数函数
给定 \(A(x)\) ，求 \(F(x)\) 使得 \(F(x)\equiv e^{A(x)}\ (mod\ x^n)\)
考虑使用牛顿迭代求解。
设 \(G(F(x))=\ln(F(x))-A(x)\) 。
但是我们还要知道 \((G(F_0(x)))'\) 的值，所以把 \(A(x)\) 当做常数，把 \(F(x)\) 当成自变量，两边求导。
现在可以套牛顿迭代了。
int exp_1[N],exp_2[N];
void Polyexp(int *A,int len){
#define sav1 exp_1
#define sav2 exp_2
int lim=1;while(lim<len) lim<<=1;
sav1[0]=1;
for(int l=2;l<=lim;l<<=1){
copy(sav2,sav1,0,l>>1),Polyln(sav2,l);
for(int i=0;i<l;++i) sav2[i]=(A[i]-sav2[i]+mo)%mo;
sav2[0]=(sav2[0]+1)%mo;
Polymul(sav1,sav2,l,l);
}
copy(A,sav1,0,len),clear(sav1,0,lim),clear(sav2,0,lim);
#undef sav1
#undef sav2
}
多项式快速幂
直接套快速幂的话会 log 方的复杂度，但我们有了指数和对数，就可以将其优化。
void Polypow(int *A,int len,int k){
Polyln(A,len);
for(int i=0;i<len;++i) A[i]=A[i]*k%mo;
Polyexp(A,len);
}
多项式开根
求 \(F(x)\) ，使得 \(F(x)\equiv \sqrt{A(x)}\ (mod\ x^n)\)
同样考虑牛顿迭代。
int sqrt_1[N],sqrt_2[N];
void Polysqrt(int *A,int len){
#define sav1 sqrt_1
#define sav2 sqrt_2
int lim=1;while(lim<len) lim<<=1;
sav1[0]=1;
for(int l=2;l<=lim;l<<=1){
for(int i=0;i<(l<<1);++i) sav2[i]=(sav1[i]<<1)%mo;
Polyinv(sav2,l);
NTT(sav1,l,1),mul(sav1,sav1,0,l),NTT(sav1,l,0);
for(int i=0;i<l;++i) sav1[i]=(sav1[i]+A[i])%mo;
Polymul(sav1,sav2,l,l);
}
copy(A,sav1,0,len),clear(sav1,0,lim),clear(sav2,0,lim);
#undef sav1
#undef sav2
}
最后来个全的！
#include<iostream>
#include<stdio.h>
#define N
#define int long long
using namespace std;
inline int read(){
int x=0,f=1;
char ch=getchar();
while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
return x*f;
}
#define mo 998244353
#define g 3
#define invg 332748118
int cir[N],rt;
int power(int x,int b){
int res=1;
while(b){
if(b&1) res=res*x%mo;
b>>=1;
x=x*x%mo;
}
return res;
}
void clear(int *A,int x,int y){
for(int i=x;i<y;++i) A[i]=0;
}
void copy(int *A,int *B,int x,int y){
for(int i=x;i<y;++i) A[i]=B[i];
}
void mul(int *A,int *B,int x,int y){
for(int i=x;i<y;++i) A[i]=A[i]*B[i]%mo;
}
void NTT(int *f,int lenth,int tag){
if(lenth!=rt)for(int i=0;i<lenth;++i)cir[i]=(cir[i>>1]>>1)|((i&1)?lenth>>1:0);
rt=lenth;
for(int i=0;i<lenth;++i) if(i<cir[i]) swap(f[i],f[cir[i]]);
for(int len=2;len<=lenth;len<<=1){
int sze=len>>1,buf=power((tag?g:invg),(mo-1)/len);
for(int i=0;i<lenth;i+=len){
int bas=1;
for(int j=i;j<i+sze;j++,bas=bas*buf%mo){
int tmp=bas*f[j+sze]%mo;
f[j+sze]=(f[j]-tmp+mo)%mo;
f[j]=(f[j]+tmp)%mo;
}
}
}
if(tag==1) return;
int invn=power(lenth,mo-2);
for(int i=0;i<lenth;++i) f[i]=f[i]*invn%mo;
}
int tmul[N];
void Polymul(int *A,int *B,int len1,int len2){
#define sav tmul
int len=1;
while(len<(len1<<1)) len<<=1;
copy(sav,B,0,len);
NTT(A,len,1),NTT(sav,len,1);
for(int i=0;i<len;++i) A[i]=A[i]*sav[i]%mo;
NTT(A,len,0);
clear(A,len2,len),clear(sav,0,len);
#undef sav
}
int inv1[N],inv2[N],inv3[N];
void Polyinv(int *A,int len){
#define sav1 inv1
#define sav2 inv2
#define sav3 inv3
int lim=1;while(lim<len) lim<<=1;
sav1[0]=power(A[0],mo-2);
for(int l=2;l<=lim;l<<=1){
int r=l<<1;copy(sav3,A,0,l);
for(int i=0;i<(l>>1);++i) sav2[i]=(sav1[i]<<1)%mo;
NTT(sav1,r,1),mul(sav1,sav1,0,r);
NTT(sav3,r,1),mul(sav1,sav3,0,r);
NTT(sav1,r,0),clear(sav1,l,r);
for(int i=0;i<l;++i) sav1[i]=(sav2[i]-sav1[i]+mo)%mo;
}
copy(A,sav1,0,len);
clear(sav1,0,lim<<1),clear(sav2,0,lim<<1),clear(sav3,0,lim<<1);
#undef sav1
#undef sav2
#undef sav3
}
int invcnt=1,inv[N];
inline void get_inv(int len){
inv[1]=1; if(invcnt>=len) return;
for(int i=invcnt+1;i<=len;++i) inv[i]=(mo-mo/i)*inv[mo%i]%mo;
}
inline void Polyder(int *A,int len){
for(int i=1;i<len;i++) A[i-1]=A[i]*i%mo; A[len-1]=0;
}
inline void Polyint(int *A,int len){
get_inv(len);
for(int i=len;i>=1;i--) A[i]=A[i-1]*inv[i]%mo; A[0]=0;
}
int ln1[N];
inline void Polyln(int *A,int len){
#define sav ln1
copy(sav,A,0,len),Polyder(A,len),Polyinv(sav,len);
Polymul(A,sav,len,len),Polyint(A,len-1),clear(sav,0,len);
#undef sav
}
int exp_1[N],exp_2[N];
void Polyexp(int *A,int len){
#define sav1 exp_1
#define sav2 exp_2
int lim=1;while(lim<len) lim<<=1;
sav1[0]=1;
for(int l=2;l<=lim;l<<=1){
copy(sav2,sav1,0,l>>1),Polyln(sav2,l);
for(int i=0;i<l;++i) sav2[i]=(A[i]-sav2[i]+mo)%mo;
sav2[0]=(sav2[0]+1)%mo;
Polymul(sav1,sav2,l,l);
}
copy(A,sav1,0,len),clear(sav1,0,lim),clear(sav2,0,lim);
#undef sav1
#undef sav2
}
void Polypow(int *A,int len,int k){
Polyln(A,len);
for(int i=0;i<len;++i) A[i]=A[i]*k%mo;
Polyexp(A,len);
}
int sqrt_1[N],sqrt_2[N];
void Polysqrt(int *A,int len){
#define sav1 sqrt_1
#define sav2 sqrt_2
int lim=1;while(lim<len) lim<<=1;
sav1[0]=1;
for(int l=2;l<=lim;l<<=1){
for(int i=0;i<(l<<1);++i) sav2[i]=(sav1[i]<<1)%mo;
Polyinv(sav2,l);
NTT(sav1,l,1),mul(sav1,sav1,0,l),NTT(sav1,l,0);
for(int i=0;i<l;++i) sav1[i]=(sav1[i]+A[i])%mo;
Polymul(sav1,sav2,l,l);
}
copy(A,sav1,0,len),clear(sav1,0,lim),clear(sav2,0,lim);
#undef sav1
#undef sav2
}
signed main(){
return 0;
}
