多项式乘法
FFT
因为有浮点数参与运算,所以可能会出现精度的问题
#include<cstdio>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const double pi=acos(-1);
const int maxn=4e6+5;
struct fs{
double x,y;
friend fs operator + (rg const fs& A,rg const fs& B){
return (fs){A.x+B.x,A.y+B.y};
}
friend fs operator - (rg const fs& A,rg const fs& B){
return (fs){A.x-B.x,A.y-B.y};
}
friend fs operator * (rg const fs& A,rg const fs& B){
return (fs){A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
}
friend fs operator / (rg const fs& A,rg const double& B){
return (fs){A.x/B,A.y/B};
}
}a[maxn],b[maxn],w[25][maxn];
int wz[maxn],n,m;
void fft(rg fs A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg fs x=A[j+k],y=A[j+k+len]*w[t0][k];
A[j+k]=x+y;
A[j+k+len]=x-y;
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
for(rg int i=0;i<lim;i++) A[i]=A[i]/(double)lim;
}
}
int main(){
n=read(),m=read();
for(rg int i=0;i<=n;i++) a[i].x=read();
for(rg int i=0;i<=m;i++) b[i].x=read();
rg int lim=1,bit=0;
for(;lim<=n+m;lim<<=1) bit++;
for(rg int i=0;i<lim;i++){
wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
}
for(rg int i=0,len=1;len<lim;len<<=1,i++){
for(rg int j=0;j<len;j++){
w[i][j]=(fs){cos(pi*j/len),sin(pi*j/len)};
}
}
fft(a,lim,1),fft(b,lim,1);
for(rg int i=0;i<lim;i++) a[i]=a[i]*b[i];
fft(a,lim,-1);
for(rg int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x+0.5));
printf("
");
return 0;
}
NTT
用原根代替复数
但必须保证模数 (p = a cdot 2^k + 1)
#include<cstdio>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e6+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,m,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int main(){
n=read(),m=read();
for(rg int i=0;i<=n;i++) a[i]=read();
for(rg int i=0;i<=m;i++) b[i]=read();
rg int lim=1,bit=0;
for(;lim<=n+m;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0,k=1;k<lim;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
ntt(a,lim,1),ntt(b,lim,1);
for(rg int i=0;i<lim;i++) a[i]=mulmod(a[i],b[i]);
ntt(a,lim,-1);
for(rg int i=0;i<=n+m;i++) printf("%d ",a[i]);
printf("
");
return 0;
}
(FFT) 和 (NTT) 预处理单位根要快不少
预处理单位根有两种写法,一种是 (O(nlogn)) 的,另一种是 (O(n)) 的
第一种因为访问内存连续实测更快
任意模数多项式乘法(MTT)
当多项式乘法进行取模的数并不是我们喜欢的模数时,就要用到 (MTT)
实现的方法有很多,这里介绍一种只需要 (2) 次 (DFT) 和 (2) 次 (IDFT) 的优秀做法
首先,我们对原来的多项式拆系数
令(k=sqrt(mod)),则
即拆成了 (A(x)/k) 和 (A(x)\%k) 两部分
那么
这样就有一个(4) 次 (DFT) 和 (4) 次 (IDFT) 的暴力做法
考虑优化,设
由于 (A, B) 的虚部都为 (0)
(P,Q) 的每一项系数都互为共轭,同样每一个点值也互为共轭
那么只需对 (P) 做一次 (DFT) ,就可以通过共轭 (O(n)) 求出 (Q) 的点值表示
具体地说,在点值表示中 (P) 的第 (k) 项与 (Q) 的第 (n-k) 项是共轭
要特判下标为 (0) 的情况
求出了 (P) 和 (Q) ,我们就能用二元一次方程的知识推出 (A_0) 和 (A_1)
对于 (B) 也是同理
此时,我们进行了两次 (DFT) 得到了所有多项式的点值表示
下面就是怎么把它们变成系数表示
仍然借鉴上面的思路,构造两个多项式
(P(x)=A_0(x)B_0(x)+iA_1(x)B_0(x))
(Q(x)=A_0(x)B_1(x)+iA_1(x)B_1(x))
通过已知的点值求出此时 (P, Q) 的点值
然后分别对 (P, Q) 做 (IDFT)
由于 (A_0 B_0,A_0 B_1,A_1 B_0,A_1 B_1)这四个多项式卷起来后的系数表示中虚部一定为 (0)
那么此时 (P) 的实部和虚部就分别为 (A_0(x) B_0(x))和 (A_1(x) B_0(x))
同样 (Q) 的实部和虚部就分别为 (A_0(x) B_1(x)) 和 (A_1(x) B_1(x))。
这样,我们又进行了两次 (IDFT) 得到了所有多项式的系数表示
代码
#include<cstdio>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5;
const double pi=acos(-1);
struct fs{
double x,y;
friend fs operator + (rg const fs& A,rg const fs& B){
return (fs){A.x+B.x,A.y+B.y};
}
friend fs operator - (rg const fs& A,rg const fs& B){
return (fs){A.x-B.x,A.y-B.y};
}
friend fs operator * (rg const fs& A,rg const fs& B){
return (fs){A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
}
friend fs operator / (rg const fs& A,rg const double& B){
return (fs){A.x/B,A.y/B};
}
}a0[maxn],b0[maxn],a1[maxn],b1[maxn],w[25][maxn],p[maxn],q[maxn];
const fs I=(fs){0,1};
int wz[maxn],n,m,mod,blo;
void fft(rg fs A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg fs x=A[j+k],y=A[j+k+len]*w[t0][k];
A[j+k]=x+y;
A[j+k+len]=x-y;
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
for(rg int i=0;i<lim;i++) A[i]=A[i]/(double)lim;
}
}
void fftfft(rg fs A[],rg fs B[],rg int lim){
for(rg int i=0;i<lim;i++){
A[i]=A[i]+I*B[i];
}
fft(A,lim,1);
B[0]=(fs){A[0].x,-A[0].y};
for(rg int i=1;i<lim;i++) B[i]=(fs){A[lim-i].x,-A[lim-i].y};
rg fs t1,t2;
for(rg int i=0;i<lim;i++){
t1=A[i],t2=B[i];
A[i]=(t1+t2)/2.0;
B[i]=(t2-t1)*I/2.0;
}
}
long long get(rg double now){
if(now>=0) return (long long)(now+0.5)%mod;
else return (long long)(now-0.5)%mod;
}
int main(){
n=read(),m=read(),mod=read();
blo=sqrt(mod)+1;
rg int aa;
for(rg int i=0;i<=n;i++){
aa=read();
aa%=mod;
a0[i].x=aa/blo;
a1[i].x=aa%blo;
}
for(rg int i=0;i<=m;i++){
aa=read();
aa%=mod;
b0[i].x=aa/blo;
b1[i].x=aa%blo;
}
rg int bit=0,lim=1;
for(;lim<=n+m;lim<<=1) bit++;
for(rg int i=0;i<lim;i++){
wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
}
for(rg int i=0,len=1;len<lim;len<<=1,i++){
for(rg int j=0;j<len;j++){
w[i][j]=(fs){cos(pi*j/len),sin(pi*j/len)};
}
}
fftfft(a0,a1,lim),fftfft(b0,b1,lim);
for(rg int i=0;i<lim;i++){
p[i]=a0[i]*b0[i]+I*a1[i]*b0[i];
q[i]=a0[i]*b1[i]+I*a1[i]*b1[i];
}
fft(p,lim,-1),fft(q,lim,-1);
for(rg int i=0;i<=n+m;i++){
long long ans=1LL*blo*blo%mod*get(p[i].x)%mod+1LL*blo*get(q[i].x+p[i].y)%mod+get(q[i].y)%mod;
printf("%lld ",ans%mod);
}
return 0;
}
多项式乘法逆
如果 (F(x)) 只剩一项
显然 (G_0) 就是 (F_0) 的逆元
如果有多项,可以递归求解
已知 (F(x)H(x) equiv 1 pmod{x^{lceil frac{n}{2} ceil}})
又因为 (F(x)G(x) equiv 1 pmod{x^{lceil frac{n}{2} ceil}})
那么 (F(x)(G(x)-H(x))equiv 1 pmod{x^{lceil frac{n}{2} ceil}})
即 (G(x)-H(x)equiv 0 pmod{x^{lceil frac{n}{2} ceil}})
两边同时平方
(G(x)^2+H(x)^2-2G(x)H(x)equiv 0 pmod{x^n})
同时乘上 (F(x))
(G(x)+F(x)H(x)^2-2H(x)equiv 0 pmod{x^n})
所以
(G(x)=2H(x)-F(x)H(x)^2(mod{x^n}))
代码
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
for(rg int i=0;i<n;i++) a[i]=read();
getinv(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("
");
return 0;
}
多项式开根
如果 (F(x)) 只剩一项
显然如果 (F_0=1), (G_0) 就是 (1)
如果有多项,可以递归求解
已知 (H(x)^2=F(x) pmod{x^{lceil frac{n}{2} ceil}})
又因为 (G(x)^2=F(x) pmod{x^{lceil frac{n}{2} ceil}})
那么 (G(x)=H(x)pmod{x^{lceil frac{n}{2} ceil}})
即 (G(x)-H(x)equiv 0 pmod{x^{lceil frac{n}{2} ceil}})
两边同时平方
(G(x)^2+H(x)^2-2G(x)H(x)equiv 0 pmod{x^n})
(F(x)+H(x)^2-2G(x)H(x)equiv 0 pmod{x^n})
所以
(G(x)=frac{F(x)+H(x)^2}{H(x)}(mod{x^n}))
套一个多项式求逆即可
代码
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
memset(B,0,maxn<<2);
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int fsqrt1[maxn],fsqrt2[maxn];
void getsqrt(rg int A[],rg int B[],rg int deg){
if(deg==1){
B[0]=1;
return;
}
getsqrt(A,B,(deg+1)>>1);
getinv(B,fsqrt1,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) fsqrt2[i]=A[i];
for(rg int i=deg;i<lim;i++) fsqrt2[i]=0;
ntt(fsqrt1,lim,1),ntt(fsqrt2,lim,1);
for(rg int i=0;i<lim;i++) fsqrt1[i]=mulmod(fsqrt1[i],fsqrt2[i]);
ntt(fsqrt1,lim,-1);
rg int ny=ksm(2,mod-2);
for(rg int i=0;i<lim;i++) B[i]=mulmod(addmod(fsqrt1[i],B[i]),ny);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
for(rg int i=0;i<n;i++) a[i]=read();
getsqrt(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("
");
return 0;
}
多项式求导
(F'(x)=frac{F(x+dx)-F(x)}{dx}=sum^{n}_{i=1}a_iinom{i}{1}x^{i-1})
因为 (dx) 趋进于无穷小,所以后面的几项都变成了 (0)
(F'(x)=sum^{n}_{i=0}a_{i + 1}(i+1)x^i)
代码
void getdao(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i-1]=mulmod(A[i],i);
}
B[deg-1]=0;
}
多项式积分
求导的逆运算
代码
void getfen(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i]=mulmod(A[i-1],ksm(i,mod-2));
}
B[0]=0;
}
多项式对数函数(多项式 ln)
对两边同时求导
(G(x)=F(A(x)),F(x)=ln(x))
(G'(x)=F'(A(x))A'(x)=frac{A'(x)}{A(x)})
所以对 (A) 求导求个逆在积一下分就算出了 (G)
代码
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
void getdao(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i-1]=mulmod(A[i],i);
}
B[deg-1]=0;
}
void getfen(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i]=mulmod(A[i-1],ksm(i,mod-2));
}
B[0]=0;
}
int fln1[maxn],fln2[maxn];
void getln(rg int A[],rg int B[],rg int deg){
memset(fln1,0,maxn<<2);
memset(fln2,0,maxn<<2);
getdao(A,fln1,deg);
getinv(A,fln2,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
ntt(fln1,lim,1),ntt(fln2,lim,1);
for(rg int i=0;i<lim;i++) fln1[i]=mulmod(fln1[i],fln2[i]);
ntt(fln1,lim,-1);
getfen(fln1,B,deg);
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
for(rg int i=0;i<n;i++) a[i]=read();
getln(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("
");
return 0;
}
多项式指数函数(多项式 exp)
(F(x)equiv F_0(x)(1-ln F_0(x)+A(x))pmod{x^n})
(F_0(x)) 指的是在模 (x^{nover 2}) 意义下的答案
代码
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
void getdao(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i-1]=mulmod(A[i],i);
}
B[deg-1]=0;
}
void getfen(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i]=mulmod(A[i-1],ksm(i,mod-2));
}
B[0]=0;
}
int fln1[maxn],fln2[maxn];
void getln(rg int A[],rg int B[],rg int deg){
memset(fln1,0,maxn<<2);
memset(fln2,0,maxn<<2);
getdao(A,fln1,deg);
getinv(A,fln2,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
ntt(fln1,lim,1),ntt(fln2,lim,1);
for(rg int i=0;i<lim;i++) fln1[i]=mulmod(fln1[i],fln2[i]);
ntt(fln1,lim,-1);
getfen(fln1,B,deg);
}
int fexp1[maxn],fexp2[maxn];
void getexp(rg int A[],rg int B[],rg int deg){
if(deg==1){
B[0]=1;
return;
}
getexp(A,B,(deg+1)>>1);
getln(B,fexp1,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) fexp2[i]=A[i];
for(rg int i=deg;i<lim;i++) fexp2[i]=0;
ntt(fexp1,lim,1),ntt(fexp2,lim,1),ntt(B,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],addmod(delmod(1,fexp1[i]),fexp2[i]));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
for(rg int i=0;i<n;i++) a[i]=read();
getexp(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("
");
return 0;
}
多项式快速幂
(G(x)=(F(x))^kpmod{x^n})
取下对数
$ln G(x)=kln F(x)pmod{x^n} $
再取一下指数
$G(x)=e^{kln F(x)}pmod{x^n} $
代码
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
void getdao(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i-1]=mulmod(A[i],i);
}
B[deg-1]=0;
}
void getfen(rg int A[],rg int B[],rg int deg){
for(rg int i=1;i<deg;i++){
B[i]=mulmod(A[i-1],ksm(i,mod-2));
}
B[0]=0;
}
int fln1[maxn],fln2[maxn];
void getln(rg int A[],rg int B[],rg int deg){
memset(fln1,0,maxn<<2);
memset(fln2,0,maxn<<2);
getdao(A,fln1,deg);
getinv(A,fln2,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
ntt(fln1,lim,1),ntt(fln2,lim,1);
for(rg int i=0;i<lim;i++) fln1[i]=mulmod(fln1[i],fln2[i]);
ntt(fln1,lim,-1);
getfen(fln1,B,deg);
}
int fexp1[maxn],fexp2[maxn];
void getexp(rg int A[],rg int B[],rg int deg){
if(deg==1){
B[0]=1;
return;
}
getexp(A,B,(deg+1)>>1);
getln(B,fexp1,deg);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) fexp2[i]=A[i];
for(rg int i=deg;i<lim;i++) fexp2[i]=0;
ntt(fexp1,lim,1),ntt(fexp2,lim,1),ntt(B,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],addmod(delmod(1,fexp1[i]),fexp2[i]));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int fmi[maxn],nowk;
void getmi(rg int A[],rg int B[],rg int deg){
getln(A,fmi,n);
for(rg int i=0;i<n;i++) fmi[i]=mulmod(fmi[i],nowk);
getexp(fmi,B,n);
}
char s[maxn];
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read();
scanf("%s",s+1);
rg int len=strlen(s+1);
for(rg int i=1;i<=len;i++){
nowk=addmod(s[i]-'0',mulmod(nowk,10));
}
for(rg int i=0;i<n;i++) a[i]=read();
getmi(a,b,n);
for(rg int i=0;i<n;i++) printf("%d ",b[i]);
printf("
");
return 0;
}
多项式除法
(F_R(x)) 等于 (F(x)) 系数翻转得到的多项式
(F(x)=Q(x)G(x)+R(x))
(F(frac{1}{x})=Q(frac{1}{x})G(frac{1}{x})+R(frac{1}{x}))
(x^n F(frac{1}{x})=x^{n-m} Q(frac{1}{x}) x^m G(frac{1}{x}) + x^{n-m+1}x^{m-1} R(frac{1}{x}))
(F_R(x)=Q_R(x)G_R(x)+x^{n-m+1}R_R(x))
(F_R(x) equiv Q_R(x)G_R(x)+x^{n-m+1}R_R(x)pmod {x^{n-m+1}})
(F_R(x) equiv Q_R(x)G_R(x)pmod {x^{n-m+1}})
(Q_R(x) equiv F_R(x)G_R^{-1}(x)pmod {x^{n-m+1}})
求一遍 (G_R) 的逆,然后就可以利用多项式乘法求出 (Q)
根据
(R(x) = F(x) - G(x)Q(x))
直接计算即可
代码
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5,mod=998244353,G=3;
int a[maxn],b[maxn],wz[maxn],n,w[23][maxn];
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
void ntt(rg int A[],rg int lim,rg int typ){
for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
for(rg int j=0,now=len<<1;j<lim;j+=now){
for(rg int k=0;k<len;k++){
rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
A[j+k]=addmod(x,y);
A[j+k+len]=delmod(x,y);
}
}
}
if(typ==-1){
std::reverse(A+1,A+lim);
rg int ny=ksm(lim,mod-2);
for(rg int i=0;i<lim;i++){
A[i]=mulmod(A[i],ny);
}
}
}
int finv[maxn];
void getinv(rg int A[],rg int B[],rg int deg){
if(deg==1){
B[0]=ksm(A[0],mod-2);
return;
}
getinv(A,B,(deg+1)>>1);
rg int lim=1,bit=0;
for(;lim<deg+deg;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg;i++) finv[i]=A[i];
for(rg int i=deg;i<lim;i++) finv[i]=0;
ntt(B,lim,1),ntt(finv,lim,1);
for(rg int i=0;i<lim;i++) B[i]=mulmod(B[i],delmod(2,mulmod(B[i],finv[i])));
ntt(B,lim,-1);
for(rg int i=deg;i<lim;i++) B[i]=0;
}
int m,fdiv1[maxn],fdiv2[maxn],ansdiv1[maxn],ansdiv2[maxn];
void getdiv(rg int A[],rg int B[],rg int deg1,rg int deg2){
std::reverse(A,A+deg1);
std::reverse(B,B+deg2);
for(rg int i=0;i<deg1-deg2+1;i++) fdiv2[i]=B[i];
getinv(fdiv2,fdiv1,deg1-deg2+1);
rg int lim=1,bit=0;
for(;lim<(deg1-deg2+1)<<1;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
for(rg int i=0;i<deg1-deg2+1;i++) fdiv2[i]=A[i];
ntt(fdiv1,lim,1),ntt(fdiv2,lim,1);
for(rg int i=0;i<lim;i++) fdiv1[i]=mulmod(fdiv1[i],fdiv2[i]);
ntt(fdiv1,lim,-1);
for(rg int i=deg1-deg2+1;i<lim;i++) fdiv1[i]=0;
memcpy(ansdiv1,fdiv1,sizeof(fdiv1));
std::reverse(ansdiv1,ansdiv1+deg1-deg2+1);
for(;lim<(deg2-1)<<1;lim<<=1) bit++;
for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
std::reverse(A,A+deg1);
std::reverse(B,B+deg2);
for(rg int i=0;i<deg2-1;i++) fdiv1[i]=ansdiv1[i],fdiv2[i]=B[i];
for(rg int i=deg2-1;i<lim;i++) fdiv1[i]=fdiv2[i]=0;
ntt(fdiv1,lim,1),ntt(fdiv2,lim,1);
for(rg int i=0;i<lim;i++) fdiv1[i]=mulmod(fdiv1[i],fdiv2[i]);
ntt(fdiv1,lim,-1);
for(rg int i=0;i<lim;i++) fdiv1[i]=delmod(A[i],fdiv1[i]);
for(rg int i=deg2-1;i<lim;i++) fdiv1[i]=0;
memcpy(ansdiv2,fdiv1,sizeof(fdiv1));
}
int main(){
for(rg int i=0,k=1;k<maxn;k<<=1,i++){
w[i][0]=1;
w[i][1]=ksm(G,(mod-1)/(k<<1));
for(rg int j=2;j<k;j++){
w[i][j]=mulmod(w[i][j-1],w[i][1]);
}
}
n=read(),m=read();
for(rg int i=0;i<=n;i++) a[i]=read();
for(rg int i=0;i<=m;i++) b[i]=read();
getdiv(a,b,n+1,m+1);
for(rg int i=0;i<n-m+1;i++) printf("%d ",ansdiv1[i]);
printf("
");
for(rg int i=0;i<m;i++) printf("%d ",ansdiv2[i]);
printf("
");
return 0;
}