多项式八分之三家桶
因为目前半家桶里只会俩所以就成四分之一了。
UPD:又会了一个,现在是八分之三了。
UUPD:虽然四个齐了,但你不感觉八分之三很带劲吗(
UUUPD:就当多项式全家桶写就完了(虽然不可能写完
多项式乘法逆
设多项式 \(A(x)\) 为 \(n\) 次多项式,求多项式 \(B\) ,使 \(A(x)*B(x)\equiv 1\pmod{x^n}\) 。
设 \(A(x)*C(x)\equiv 1\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\) ,那么
递归求解。由主定理复杂度 \(O(n\log n)\) 。
洛谷P4238[模板]多项式乘法逆
##include<bits/stdc++.h>
using namespace std;
namespace IO{
typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp,int len=0){
if(x<0) putchar('-'), x=-x;
do{ outp[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
void ckmin(int& x,int y){ x=x<y?x:y; }
void ckmax(int& x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=400010,mod=998244353;
int n;
namespace Poly_Calculation{
const LL rt=3,two=998244355;
int l,len,r[NN];
LL a[NN],b[NN],g[NN],tp[NN];
void pls(LL& x,LL y){ (x+=y)>=mod?(x-=mod):(x<0?(x+=mod):0); }
LL qpow(LL a,int b=mod-2,LL res=1){
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
void prework(int n){
for(len=1,l=0;len<=n+n;len<<=1) ++l;
for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
g[0]=1; g[1]=qpow(rt,(mod-1)/len);
for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
}
void ntt(LL* a,int n){
for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
LL tmp=g[j*t]*a[i+j+d]%mod;
pls(a[i+j+d]=a[i+j],-tmp); pls(a[i+j],tmp);
}
}
void polyinv(LL* a,LL* res,int deg){
if(deg==1) return void(res[0]=qpow(a[0]));
polyinv(a,res,deg+1>>1); prework(deg);
for(int i=0;i<deg;i++) tp[i]=a[i];
for(int i=deg;i<len;i++) tp[i]=0;
ntt(res,len); ntt(tp,len);
for(int i=0;i<len;i++) res[i]=(two-res[i]*tp[i]%mod)*res[i]%mod, g[i]=qpow(g[i]);
ntt(res,len); LL inv=qpow(len);
for(int i=0;i<deg;i++) res[i]=res[i]*inv%mod;
for(int i=deg;i<len;i++) res[i]=0;
}
} using namespace Poly_Calculation;
signed main(){
n=read();
for(int i=0;i<n;i++) a[i]=read();
polyinv(a,b,n);
for(int i=0;i<n;i++) write(b[i],' ');
return puts(""),0;
}
多项式开根
已知 \(n\) 多项式 \(A(x)\) ,求多项式 \(B(x)\) ,满足 \(B^2(x)\equiv A(x)\pmod{x^n}\) 。
类似地,设 \(C^2(x)\equiv A(x)\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\) ,那么
递归并求逆求解。
这份板子因为某种神秘力量常数过丑了。有时间的话应该会来改改吧。(
UPD:取模优化去掉开O2跑快了三秒。。不懂什么原理。虽然还是挺慢但至少能看了。先这样吧
洛谷P5205[模板]多项式开根
##include<bits/stdc++.h>
using namespace std;
namespace IO{
typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp,int len=0){
if(x<0) putchar('-'), x=-x;
do{ outp[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
void ckmin(int& x,int y){ x=x<y?x:y; }
void ckmax(int& x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=400010,mod=998244353;
int n;
namespace Poly_Calculation{
const LL rt=3,two=998244355,inv2=499122177;
int l,len,r[NN];
LL a[NN],b[NN],g[NN],iv[NN],tp[NN];
void pls(LL& x,LL y){ (x+=y)>=mod?(x-=mod):(x<0?(x+=mod):0); }
LL qpow(LL a,int b=mod-2,LL res=1){
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
void prework(int n){
for(len=1,l=0;len<=(n<<1);len<<=1) ++l;
for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
g[0]=1; g[1]=qpow(rt,(mod-1)/len);
for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
}
void ntt(LL* a,int n){
for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
LL tmp=g[j*t]*a[i+j+d]%mod;
a[i+j+d]=(a[i+j]-tmp+mod)%mod;
a[i+j]=(a[i+j]+tmp)%mod;
}
}
void polyinv(LL* a,LL* res,int deg){
if(deg==1) return void(res[0]=qpow(a[0]));
polyinv(a,res,deg+1>>1); prework(deg);
for(int i=0;i<deg;i++) tp[i]=a[i];
for(int i=deg;i<len;i++) tp[i]=0;
ntt(res,len); ntt(tp,len);
for(int i=0;i<len;i++) res[i]=(two-res[i]*tp[i]%mod)*res[i]%mod, g[i]=qpow(g[i]);
ntt(res,len); LL inv=qpow(len);
for(int i=0;i<deg;i++) res[i]=res[i]*inv%mod;
for(int i=deg;i<len;i++) res[i]=0;
}
void polysqr(LL* a,LL *res,int deg){
if(deg==1) return void(res[0]=1);
polysqr(a,res,deg+1>>1);
for(int i=0;i<len;i++) iv[i]=0;
polyinv(res,iv,deg); prework(deg);
for(int i=0;i<deg;i++) tp[i]=a[i];
for(int i=deg;i<len;i++) tp[i]=0;
ntt(res,len); ntt(iv,len); ntt(tp,len);
for(int i=0;i<len;i++) res[i]=(tp[i]*iv[i]%mod+res[i])*inv2%mod, g[i]=qpow(g[i]);
ntt(res,len); LL inv=qpow(len);
for(int i=0;i<deg;i++) res[i]=res[i]*inv%mod;
for(int i=deg;i<len;i++) res[i]=0;
}
} using namespace Poly_Calculation;
signed main(){
n=read();
for(int i=0;i<n;i++) a[i]=read();
polysqr(a,b,n);
for(int i=0;i<n;i++) write(b[i],' ');
return puts(""),0;
}
多项式对数函数
已知 \(n\) 次多项式 \(A\) ,求多项式 \(B\) ,满足 \(B(x)\equiv \ln(A(x))\pmod{x^n}\) 。
\(\ln\) 很ex,求导化掉。
求乘法逆之后求出导数再积分回去。
不定积分和求导是互逆的,
洛谷P4725[模板]多项式对数函数
##include<bits/stdc++.h>
using namespace std;
namespace IO{
typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp,int len=0){
if(x<0) putchar('-'), x=-x;
do{ outp[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
void ckmin(int& x,int y){ x=x<y?x:y; }
void ckmax(int& x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=600010,mod=998244353;
int n;
namespace Poly_Calculation{
const LL rt=3,two=mod+2;
int l,len,r[NN];
LL a[NN],b[NN],g[NN],iv[NN],tp[NN];
LL qpow(LL a,int b=mod-2,LL res=1){
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
void prework(int n){
for(len=1,l=0;len<=(n<<1);len<<=1) ++l;
for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
g[0]=1; g[1]=qpow(rt,(mod-1)/len);
for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
}
void ntt(LL* a,int n){
for(int i=0;i<n;i++) if(r[i]>i) swap(a[i],a[r[i]]);
for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
LL tmp=g[t*j]*a[i+j+d]%mod;
a[i+j+d]=(a[i+j]-tmp+mod)%mod; a[i+j]=(a[i+j]+tmp)%mod;
}
}
void deriv(LL* a,LL* res,int n){ for(int i=1;i<n;i++) res[i-1]=a[i]*i%mod; res[n-1]=0; }
void integ(LL* a,LL* res,int n){ for(int i=1;i<n;i++) res[i]=a[i-1]*qpow(i)%mod; res[0]=0; }
void polyinv(LL* a,LL* res,int deg){
if(deg==1) return void(res[0]=qpow(a[0]));
polyinv(a,res,deg+1>>1); prework(deg);
for(int i=0;i<deg;i++) tp[i]=a[i];
for(int i=deg;i<len;i++) tp[i]=0;
ntt(tp,len); ntt(res,len);
for(int i=0;i<len;i++) res[i]=(two-res[i]*tp[i]%mod)*res[i]%mod, g[i]=qpow(g[i]);
ntt(res,len); LL inv=qpow(len);
for(int i=0;i<deg;i++) res[i]=res[i]*inv%mod;
for(int i=deg;i<len;i++) res[i]=0;
}
void polyln(LL* a,LL* res,int deg){
polyinv(a,iv,deg); prework(deg);
for(int i=0;i<len;i++) tp[i]=0;
deriv(a,tp,deg);
ntt(iv,len); ntt(tp,len);
for(int i=0;i<len;i++) tp[i]=tp[i]*iv[i]%mod, g[i]=qpow(g[i]);
ntt(tp,len); LL inv=qpow(len);
for(int i=0;i<len;i++) tp[i]=tp[i]*inv%mod;
integ(tp,res,deg);
}
} using namespace Poly_Calculation;
signed main(){
n=read();
for(int i=0;i<n;i++) a[i]=read();
polyln(a,b,n);
for(int i=0;i<n;i++) write(b[i],' ');
return puts(""),0;
}
多项式指数函数
已知 \(n\) 次多项式 \(A(x)\) ,求 \(n\) 次多项式 \(B(x)\) ,满足 \(B(x)\equiv e^{A(x)}\pmod{x^n}\)
仍考虑递归求解。如果已知函数 \(G(x)\) ,求满足 \(G(F(z))\equiv 0(\mod z^n)\) 的 \(n\) 次多项式 \(F(z)\) ,可以用多项式牛顿迭代。
(洛谷题解第一篇的图
现在令 \(G(x)=\ln B(x)-A(x)\) ,令 \(B_0(x)\) 为 \(\mod x^{\left\lceil\frac{n}{2}\right\rfloor}\) 意义下的解,那么就是要求 \(G\) 的一个零点。
对 \(G\) 求导,把 \(A(x)\) 看作常数,那么 \(G'(B(x))=\frac{1}{B(x)}\) ,代入多项式牛顿迭代,有
求 \(\ln\) 并递归求解。
洛谷P4726[模板]多项式指数函数
##include<bits/stdc++.h>
using namespace std;
namespace IO{
typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp,int len=0){
if(x<0) putchar('-'), x=-x;
do{ outp[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
void ckmin(int& x,int y){ x=x<y?x:y; }
void ckmax(int& x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=600010,mod=998244353;
int n;
namespace Poly_Calculation{
const LL rt=3,one=mod+1,two=mod+2;
int l,len,r[NN];
LL a[NN],g[NN],b[NN],iv[NN],tp[NN],ln[NN];
LL qpow(LL a,int b=mod-2,LL res=1){
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
void prework(int n){
for(len=1,l=0;len<(n<<1);len<<=1) ++l;
for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
g[0]=1; g[1]=qpow(rt,(mod-1)/len);
for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
}
void ntt(LL* a,int n){
for(int i=0;i<n;i++) if(r[i]>i) swap(a[i],a[r[i]]);
for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
LL tmp=g[t*j]*a[i+j+d]%mod;
a[i+j+d]=(a[i+j]-tmp+mod)%mod; a[i+j]=(a[i+j]+tmp)%mod;
}
}
void polyinv(LL* a,LL* res,int deg){
if(deg==1) return void(res[0]=qpow(a[0]));
polyinv(a,res,deg+1>>1); prework(deg);
for(int i=0;i<deg;i++) tp[i]=a[i];
ntt(tp,len); ntt(res,len);
for(int i=0;i<len;i++) res[i]=(two-res[i]*tp[i]%mod)*res[i]%mod, g[i]=qpow(g[i]);
ntt(res,len); LL inv=qpow(len);
for(int i=0;i<deg;i++) res[i]=res[i]*inv%mod, tp[i]=0;
for(int i=deg;i<len;i++) res[i]=tp[i]=0;
}
void deriv(LL* a,LL* res,int n){ for(int i=1;i<n;i++) res[i-1]=a[i]*i%mod; res[n-1]=0; }
void integ(LL* a,LL* res,int n){ for(int i=1;i<n;i++) res[i]=a[i-1]*qpow(i)%mod; res[0]=0; }
void polyln(LL* a,LL* res,int deg){
polyinv(a,iv,deg); deriv(a,tp,deg); prework(deg);
ntt(iv,len); ntt(tp,len);
for(int i=0;i<len;i++) tp[i]=tp[i]*iv[i]%mod, g[i]=qpow(g[i]);
ntt(tp,len); LL inv=qpow(len);
for(int i=0;i<deg;i++) tp[i]=tp[i]*inv%mod;
integ(tp,res,deg);
for(int i=0;i<len;i++) tp[i]=iv[i]=0;
}
void polyexp(LL* a,LL* res,int deg){
if(deg==1) return void(res[0]=1);
polyexp(a,res,deg+1>>1); polyln(res,ln,deg); prework(deg);
for(int i=0;i<deg;i++) tp[i]=a[i];
ntt(tp,len); ntt(res,len); ntt(ln,len);
for(int i=0;i<len;i++) res[i]=(one-ln[i]+tp[i])*res[i]%mod, g[i]=qpow(g[i]);
ntt(res,len); LL inv=qpow(len);
for(int i=0;i<deg;i++) res[i]=res[i]*inv%mod, tp[i]=ln[i]=0;
for(int i=deg;i<len;i++) res[i]=tp[i]=ln[i]=0;
}
} using namespace Poly_Calculation;
signed main(){
n=read();
for(int i=0;i<n;i++) a[i]=read();
polyexp(a,b,n);
for(int i=0;i<n;i++) write(b[i],' ');
return puts(""),0;
}
多项式半家桶
用 namespace
套起来能算封装吗(
洛谷P4389付公主的背包
#include<bits/stdc++.h>
#define int long long
using namespace std;
namespace IO{
typedef unsigned long long ULL;
typedef double DB; typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp,int len=0){
if(x<0) putchar('-'), x=-x;
do{ outp[len++]=(x%10)^48; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
void ckmin(int& x,int y){ x=x<y?x:y; }
void ckmax(int& x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=600010,mod=998244353;
int n,m,buc[NN],inv[NN];
int qpow(int a,int b=mod-2,int res=1){
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
namespace PolyCalculation{
const int rt=3,one=mod+1,two=mod+2,inv2=mod+1>>1;
int l,len,r[NN],a[NN],g[NN],tp[NN],iv[NN],ln[NN],res[NN];
void revg(){ g[1]=qpow(g[1]); for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod; }
void deriv(int* a,int* r,int n){ for(int i=1;i<n;i++) r[i-1]=i*a[i]%mod; r[n-1]=0; }
void integ(int* a,int* r,int n){ for(int i=1;i<n;i++) r[i]=a[i-1]*qpow(i)%mod; r[0]=0; }
void prework(int n){
for(len=1,l=0;len<=n;len<<=1) ++l;
for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
g[0]=1; g[1]=qpow(rt,(mod-1)/len);
for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
}
void ntt(int* a,int n){
for(int i=0;i<n;i++) if(r[i]>i) swap(a[i],a[r[i]]);
for(int t=n>>1,d=1;d<n;t>>=1,d<<=1)
for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
int tmp=g[t*j]*a[i+j+d]%mod;
a[i+j+d]=(a[i+j]-tmp+mod)%mod; a[i+j]=(a[i+j]+tmp)%mod;
}
}
void polyinv(int* a,int* r,int deg){
if(deg==1) return void(r[0]=qpow(a[0]));
polyinv(a,r,deg+1>>1); prework(deg<<1);
for(int i=0;i<deg;i++) tp[i]=a[i];
ntt(tp,len); ntt(r,len);
for(int i=0;i<len;i++) r[i]=(two-r[i]*tp[i]%mod)*r[i]%mod;
revg(); ntt(r,len); int inv=qpow(len);
for(int i=0;i<deg;i++) r[i]=r[i]*inv%mod, tp[i]=0;
for(int i=deg;i<len;i++) r[i]=tp[i]=0;
}
void polysqr(LL* a,LL *r,int deg){
if(deg==1) return void(r[0]=1);
polysqr(a,r,deg+1>>1); polyinv(r,iv,deg); prework(deg<<1);
for(int i=0;i<deg;i++) tp[i]=a[i];
ntt(r,len); ntt(iv,len); ntt(tp,len);
for(int i=0;i<len;i++) r[i]=(tp[i]*iv[i]%mod+r[i])*inv2%mod;
revg(); ntt(r,len); LL inv=qpow(len);
for(int i=0;i<deg;i++) r[i]=r[i]*inv%mod, iv[i]=tp[i]=0;
for(int i=deg;i<len;i++) r[i]=0, iv[i]=tp[i]=0;
}
void polyln(int* a,int* r,int deg){
polyinv(a,iv,deg); deriv(a,tp,deg);
prework(deg<<1); ntt(tp,len); ntt(iv,len);
for(int i=0;i<len;i++) tp[i]=tp[i]*iv[i]%mod;
revg(); ntt(tp,len); int inv=qpow(len);
for(int i=0;i<len;i++) tp[i]=tp[i]*inv%mod;
integ(tp,r,deg);
for(int i=0;i<len;i++) tp[i]=iv[i]=0;
}
void polyexp(int* a,int* r,int deg){
if(deg==1) return void(r[0]=1);
polyexp(a,r,deg+1>>1); polyln(r,ln,deg); prework(deg<<1);
for(int i=0;i<deg;i++) tp[i]=a[i];
ntt(tp,len); ntt(ln,len); ntt(r,len);
for(int i=0;i<len;i++) r[i]=(one-ln[i]+tp[i])*r[i]%mod;
revg(); ntt(r,len); int inv=qpow(len);
for(int i=0;i<deg;i++) r[i]=r[i]*inv%mod, ln[i]=tp[i]=0;
for(int i=deg;i<len;i++) r[i]=ln[i]=tp[i]=0;
}
} using namespace PolyCalculation;
signed main(){
n=read(); m=read(); inv[0]=inv[1]=1;
for(int i=1;i<=n;i++) ++buc[read()];
for(int i=2;i<=m;i++) inv[i]=mod-mod/i*inv[mod%i]%mod;
for(int i=1;i<=m;i++) if(buc[i])
for(int j=i,k=1;j<=m;j+=i,k++)
a[j]=(a[j]+buc[i]*inv[k]);
polyexp(a,res,m+1);
for(int i=1;i<=m;i++) write(res[i],'\n');
return 0;
}
多项式另八分之五家桶(
任意模数多项式乘法
当出题人毒瘤模数不符合 \(\tt{NTT}\) 要求时就要写这个。有两种写法,一是三模数 \(\tt{NTT}\) ,一是拆系数 \(\tt{FFT}\) 。前者要做变换的次数较多,而后者经过一系列NB优化后一共只用四次 \(\tt{FFT}\) ,因此学了后者。
三模数 \(\tt{NTT}\) 思路大体是用三个 \(\tt{NTT}\) 模数通过中国剩余定理合并出一个大于运算中最大值的超大模数,合并答案后取模。好像要 \(9\) 次 \(\tt{NTT}\) 。
拆系数 \(\tt{FFT}\) (\(\tt{MTT}\))将多项式 \(F(x)=\sum_{i=0}^na_ix^i\) 拆为 \(F(x)=\sum_{i=0}^n(b_i+c_i)x^i\) ,其中 \(b_i=\left\lfloor\frac{a_i}{\sqrt{mod}}\right\rfloor\) , \(c_i=a_i\mod\sqrt{mod}\) 。实际运算中 \(\sqrt{mod}\) 可以直接取 \(2^{15}\) 。化简式子后可以 \(7\) 次 \(\tt{FFT}\) ,再发掘一些共轭与 \(\tt{DFT}\) 的性质后可以优化到 \(4\) 次 \(\tt{FFT}\) 。具体可见洛谷题解
洛谷P4245[模板]任意模数多项式乘法
#include<bits/stdc++.h>
using namespace std;
namespace IO{
typedef long double DB;
typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp,int len=0){
if(x<0) putchar('-'), x=-x;
do{ outp[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
void ckmin(int& x,int y){ x=x<y?x:y; }
void ckmax(int& x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=800010;
int n,m,mod,f[NN],g[NN],res[NN];
namespace Poly_Calculation{
const DB PI=acos(-1.0);
struct compx{
DB r,i;
compx(){} compx(DB x,DB y){ r=x; i=y; }
compx conj(){ return compx(r,-i); }
compx operator+(const compx& t)const{ return compx(r+t.r,i+t.i); }
compx operator-(const compx& t)const{ return compx(r-t.r,i-t.i); }
compx operator*(const compx& t)const{ return compx(r*t.r-i*t.i,r*t.i+i*t.r); }
}a[NN],b[NN],c[NN],d[NN],w[NN];
int l,len,r[NN];
void prework(int n){
for(len=1,l=0;len<n;len<<=1) ++l;
for(int i=0;i<len;i++){
r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
w[i]=compx(cos(2.0*i*PI/len),sin(2.0*i*PI/len));
}
}
void fft(compx* a,int n){
for(int i=0;i<n;i++) if(r[i]>i) swap(a[i],a[r[i]]);
for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
compx tmp=w[j*t]*a[i+j+d];
a[i+j+d]=a[i+j]-tmp; a[i+j]=a[i+j]+tmp;
}
}
void mtt(int* f,int* g,int* res,int n){
prework(n);
for(int i=0;i<n;i++){
a[i]=compx(f[i]>>15,f[i]&32767);
c[i]=compx(g[i]>>15,g[i]&32767);
}
fft(a,len); fft(c,len);
for(int i=1;i<len;i++) b[i]=a[len-i].conj(), d[i]=c[len-i].conj();
b[0]=a[0].conj(); d[0]=c[0].conj();
for(int i=0;i<len;i++){
compx aa=(a[i]+b[i])*compx(0.5,0), bb=(a[i]-b[i])*compx(0,-0.5);
compx cc=(c[i]+d[i])*compx(0.5,0), dd=(c[i]-d[i])*compx(0,-0.5);
a[i]=aa*cc+compx(0,1)*(aa*dd+cc*bb); b[i]=bb*dd; w[i].i=-w[i].i;
}
fft(a,len); fft(b,len);
for(int i=0;i<len;i++){
LL x=(LL)(a[i].r/len+0.5)%mod,y=(LL)(a[i].i/len+0.5)%mod,z=(LL)(b[i].r/len+0.5)%mod;
res[i]=((1ll<<30)*x%mod+(1ll<<15)*y%mod+z+mod)%mod;
}
}
} using namespace Poly_Calculation;
signed main(){
n=read()+1; m=read()+1; mod=read();
for(int i=0;i<n;i++) f[i]=read();
for(int i=0;i<m;i++) g[i]=read();
mtt(f,g,res,n+m);
for(int i=0;i<n+m-1;i++) write(res[i],' ');
return puts(""),0;
}
分治 \(\tt{FFT}\)
求 \(f_i=\sum f_jg_{i-j}\) 这样像卷积但又有些ex的式子。做CDQ分治就好。
洛谷P4721[模板]分治FFT
#include<bits/stdc++.h>
using namespace std;
namespace IO{
#define int LL
typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp,int len=0){
if(x<0) putchar('-'), x=-x;
do{ outp[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
void ckmin(int& x,int y){ x=x<y?x:y; }
void ckmax(int& x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=600010,mod=998244353;
int n;
int qpow(int a,int b=mod-2,int res=1){
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
namespace Poly_Calculation{
const int rt=3;
int l,len,r[NN],a[NN],b[NN],g[NN],tp[NN],pt[NN];
void prework(int n){
for(len=1,l=0;len<=n;len<<=1) ++l;
for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
g[0]=1; g[1]=qpow(rt,(mod-1)/len);
for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
}
void ntt(int* a,int n){
for(int i=0;i<n;i++) if(r[i]>i) swap(a[i],a[r[i]]);
for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
int tmp=g[j*t]*a[i+j+d]%mod;
a[i+j+d]=(a[i+j]-tmp+mod)%mod; a[i+j]=(a[i+j]+tmp)%mod;
}
}
void conquer(int* ff,int* gg,int l,int r){
if(l==r){ if(!l) ff[0]=1; return; }
int mid=l+r>>1; conquer(ff,gg,l,mid); prework(r-l+1);
for(int i=l;i<=mid;i++) tp[i-l]=ff[i];
for(int i=0;i<=r-l;i++) pt[i]=gg[i];
ntt(tp,len); ntt(pt,len);
for(int i=0;i<len;i++) tp[i]=tp[i]*pt[i]%mod, g[i]=qpow(g[i]);
ntt(tp,len); int inv=qpow(len);
for(int i=mid+1;i<=r;i++) ff[i]=(ff[i]+tp[i-l]*inv)%mod;
for(int i=0;i<len;i++) tp[i]=pt[i]=0;
conquer(ff,gg,mid+1,r);
}
} using namespace Poly_Calculation;
signed main(){
n=read();
for(int i=1;i<n;i++) b[i]=read();
conquer(a,b,0,n-1);
for(int i=0;i<n;i++) write(a[i],' ');
return puts(""),0;
}
多项式除法/多项式取模
已知 \(n\) 次多项式 \(F\) 与 \(m\) 次多项式 \(G\) ,求 \(n-m\) 次多项式 \(Q\) 与次数不大于 \(m\) 的多项式 \(R\) ,满足 \(F(x)=G(x)Q(x)+R(x)\) 。
把多项式系数翻转,然后可以提出 \(x\) 的次幂,令 \(F\) 反转后为 \(F_r\) ,有
在模 \(x^{n-m+1}\) 意义下后面的东西为 \(0\) ,于是模 \(x^{n-m+1}\) 求逆后容易得出两个多项式。
洛谷P4512[模板]多项式除法
#include<bits/stdc++.h>
using namespace std;
namespace IO{
#define int LL
typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp,int len=0){
if(x<0) putchar('-'), x=-x;
do{ outp[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
void ckmin(int& x,int y){ x=x<y?x:y; }
void ckmax(int& x,int y){ x=x>y?x:y; }
} using namespace IO;
const int NN=400010,mod=998244353;
int n,m,f[NN],q[NN],res[NN];
int qpow(int a,int b=mod-2,int res=1){
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
namespace Poly_Calculation{
const int rt=3,two=mod+2;
int l,len,r[NN],g[NN],ff[NN],qq[NN],hh[NN],dd[NN],vv[NN],tp[NN],pt[NN];
void revg(){ g[1]=qpow(g[1]); for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod; }
void prework(int n){
for(len=1,l=0;len<=n;len<<=1) ++l;
for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
g[0]=1; g[1]=qpow(rt,(mod-1)/len);
for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
}
void ntt(int* a,int n){
for(int i=0;i<n;i++) if(r[i]>i) swap(a[i],a[r[i]]);
for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
int tmp=g[j*t]*a[i+j+d]%mod;
a[i+j+d]=(a[i+j]-tmp+mod)%mod; a[i+j]=(a[i+j]+tmp)%mod;
}
}
void mul(int* a,int* b,int* c,int n){
for(int i=0;i<n;i++) tp[i]=a[i], pt[i]=b[i];
prework(n); ntt(tp,len); ntt(pt,len);
for(int i=0;i<len;i++) tp[i]=tp[i]*pt[i]%mod;
revg(); ntt(tp,len); int inv=qpow(len);
for(int i=0;i<len;i++) c[i]=tp[i]*inv%mod, tp[i]=pt[i]=0;
}
void polyinv(int* a,int* res,int deg){
if(deg==1) return void(res[0]=qpow(a[0]));
polyinv(a,res,deg+1>>1); prework(deg<<1);
for(int i=0;i<deg;i++) tp[i]=a[i];
ntt(tp,len); ntt(res,len);
for(int i=0;i<len;i++) res[i]=(two-res[i]*tp[i]%mod)*res[i]%mod;
revg(); ntt(res,len); int inv=qpow(len);
for(int i=0;i<deg;i++) res[i]=res[i]*inv%mod, tp[i]=0;
for(int i=deg;i<len;i++) res[i]=tp[i]=0;
}
void polymod(int* a,int* b,int* q,int* r,int n,int m){
for(int i=0;i<=n;i++) q[n-i]=ff[i]=a[i];
for(int i=0;i<=m;i++) dd[m-i]=hh[i]=b[i];
for(int i=n-m+2;i<=m;i++) dd[i]=0;
polyinv(dd,vv,n-m+1); mul(q,vv,q,n<<1); reverse(q,q+n-m+1);
for(int i=n-m+1;i<=n;i++) q[i]=0; mul(hh,q,hh,n<<1);
for(int i=0;i<m;i++) r[i]=(ff[i]+mod-hh[i])%mod;
for(int i=0;i<len;i++) dd[i]=ff[i]=hh[i]=0;
}
} using namespace Poly_Calculation;
signed main(){
n=read(); m=read();
for(int i=0;i<=n;i++) f[i]=read();
for(int i=0;i<=m;i++) q[i]=read();
polymod(f,q,qq,res,n,m);
for(int i=0;i<=n-m;i++)write(qq[i],' ');puts("");
for(int i=0;i<m;i++) write(res[i],' '); puts("");
return 0;
}