fft && ntt 模板
作者:互联网
FFT模板
点击查看代码
#include<bits/stdc++.h>
using namespace std;
template<typename T>
inline void read(T &x){
x=0;T fl=1;char tmp=getchar();
while(tmp<'0'||tmp>'9')fl=tmp=='-'?-fl:fl,tmp=getchar();
while(tmp>='0'&&tmp<='9')x=(x<<1)+(x<<3)+tmp-'0',tmp=getchar();
x=x*fl;
}
const double Pi=acos(-1);
const int maxn=4e6;
struct FastFourierTransform{
complex<double>omega[maxn],iomega[maxn];
void init(const int n){
for(int i=0;i<n;i++){
omega[i]=complex<double>(cos(2*Pi/n*i),sin(2*Pi/n*i));
iomega[i]=conj(omega[i]);
}
}
void transform(complex<double>*a,int n,complex<double> *omega){
int k=0;
while((1<<k)<n)k++;
for(int i=0;i<n;i++){
int t=0;
for(int j=0;j<k;j++)
if(i&(1<<j))t|=1<<k-j-1;
if(t<i)swap(a[i],a[t]);
}
for(int l=2;l<=n;l<<=1){
int m=l/2;
for(complex<double> *p=a;p!=a+n;p+=l){
for(int i=0;i<m;i++){
complex<double> t=omega[n/l*i]*p[i+m];
p[i+m]=p[i]-t;
p[i]=p[i]+t;
}
}
}
}
void dft(complex<double>*a,const int n){
transform(a,n,omega);
}
void idft(complex<double>*a,const int n){
transform(a,n,iomega);
for(int i=0;i<n;i++)
a[i]/=n;
}
}fft;
inline int multiply(int *a1,int n1,int *a2,const int n2){
int n=1;
while(n<n1+n2)n<<=1;
complex<double>c1[maxn],c2[maxn];
for(int i=0;i<n1;i++)
c1[i]=a1[i];
for(int i=0;i<n2;i++)
c2[i]=a2[i];
fft.init(n);
fft.dft(c1,n),fft.dft(c2,n);
for(int i=0;i<n;i++) c1[i]=c1[i]*c2[i];
fft.idft(c1,n);
n1=n1+n2-1;
for(int i=0;i<n1+n2-1;i++) a1[i]=floor(c1[i].real()+0.5);
return n1;
}
int n,m;
int a[maxn],b[maxn];
signed main(){
cin>>n>>m;n++,m++;
for(int i=0;i<n;i++)
read(a[i]);
for(int i=0;i<m;i++)
read(b[i]);
n=multiply(a,n,b,m);
for(int i=0;i<n;i++)
printf("%d ",a[i]);
return 0;
}
NTT模板
求原根,可以先求出p-1的质因子。由于原根往往不大,可以从2到p-1依次枚举i,
\(若i^{(p-1)/p_i}≠1(mod p),{\forall}p_i为(p-1)的质因数,则i为p的一个原根。\)
如此将复数单位根替换成模p意义下的
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=998244353;
const int maxn=2.5e6;
template<typename T>
inline void read(T &x){
x=0;T fl=1;char tmp=getchar();
while(tmp<'0'||tmp>'9')fl=tmp=='-'?-fl:fl,tmp=getchar();
while(tmp>='0'&&tmp<='9')x=(x<<1)+(x<<3)+tmp-'0',tmp=getchar();
x=x*fl;
}
inline ll pw(ll x,ll n,ll p){
ll ans=1;
while(n){
if(n&1)ans=ans*x%p;
x=x*x%p,n>>=1;
}
return ans;
}
inline ll root(const ll p){
ll pri[60],cnt=0;
ll x=p-1;
for(int k=2;k*k<=p-1;k++){
if(x%k==0){
pri[++cnt]=k;
while(x%k==0)x/=k;
}
}
if(x>1)pri[++cnt]=x;
int fl;
for(int i=2;i<=p;i++){
fl=0;
for(int j=1;j<=cnt;j++){
if(pw(i,(p-1)/pri[j],p)==1){
fl=1;
break;
}
}
if(!fl)return i;
}
throw;
}
inline void exgcd(const ll a,const ll b,ll &x, ll &y){
if(!b)x=1,y=0;
else exgcd(b,a%b,y,x),y-=a/b*x;
}
inline ll inv(const ll a,const ll p){//must exist
ll x,y;
exgcd(a,p,x,y);
return (x%p+p)%p;
}
struct NumberTheoreticTransform{
ll omega[maxn],iomega[maxn];
void init(const int n){
ll g=root(mod),x=pw(g,(mod-1)/n,mod);
omega[0]=iomega[0]=1;
for(int i=1;i<n;i++){
omega[i]=omega[i-1]*x;
iomega[i]=inv(omega[i],mod);
}
}
void transform(ll *a,const int n,ll *omega){
int k=0;
while((1<<k)<n)k++;
for(int i=0;i<n;i++){
int t=0;
for(int j=0;j<k;j++)if(i&(1<<j))t|=1<<k-j-1;
if(i>t)swap(a[i],a[t]);
}
for(int l=2;l<=n;l++){
int m=l/2;
for(ll *p=a;p!=a+n;p+=l)
for(int i=0;i<m;i++){
ll t=omega[n/l*i]*p[i+m]%mod;
p[i+m]=(p[i]-t+mod)%mod;
p[i]=(p[i]+t)%mod;
}
}
}
void dft(ll *a,const int n){
transform(a,n,omega);
}
void idft(ll *a,const int n){
transform(a,n,iomega);
ll x=inv(n,mod);
for(int i=0;i<n;i++)
a[i]=a[i]*x%mod;
}
}ntt;
标签:tmp,const,int,ll,fft,ntt,&&,fl,void 来源: https://www.cnblogs.com/xyc1719/p/16533688.html