# 正紧地多项式求exp 牛客练习赛24F题解

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define re register
#define rep(i,s,t) for(re int i=s;i<=t;++i)
#define _rep(i,s,t) for(re int i=s;i>=t;--i)
#define Rep(i,s,t) for(re int i=s;i<t;++i)
#define go(x) for(re int e=las[x];e;e=nxt[e])
#define re register
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define pii pair<int,int>
#define pi acos(-1)
#define ms(f,x) memset(f,x,sizeof f)
#define open(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
namespace IO{
#define gc getchar()
#define pc(x) putchar(x)
template<typename T>inline void read(T &x){
x=0;int f=1;char ch=gc;while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=gc;}
while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+ch-'0',ch=gc;x*=f;return;
}
template<typename T>inline void write(T x=0){
T wr[51];wr[0]=0;if(x<0)pc('-'),x=-x;if(!x)pc(48);
while(x)wr[++wr[0]]=x%10,x/=10;while(wr[0])pc(48+wr[wr[0]--]);return;
}
}
using IO::write;
using namespace std;
typedef long long ll;
const int N=1e6+11,mod=19260817;
int n,m;
int a[N],b[N],p[N],c[N],
d[N],e[N],f[N],inv[N];
int lg2;
typedef double db;
struct cp{
db r,i;
cp(db a=0.,db b=0.){r=a,i=b;}
inline cp operator +(cp A)const{return cp(r+A.r,i+A.i);}
inline cp operator -(cp A)const{return cp(r-A.r,i-A.i);}
inline cp operator *(cp A)const{return cp(r*A.r-i*A.i,r*A.i+i*A.r);}
inline cp operator *(int A)const{return cp(A,0);}
}a1[N],a2[N],b1[N],b2[N],c1[N],c2[N],c3[N];
inline int inc(int x,int y){
re int res=x+y;
if(res>=mod)res-=mod;
if(res<0)res+=mod;
return res;
}
inline int fp(int a,int b){
if(b>=mod-1)b-=mod-1;
if(b<0)b+=mod-1;
re int res=1;
for(;b;b>>=1,a=1ll*a*a%mod)
if(b&1)
res=1ll*res*a%mod;
return res;
}
inline void fft(cp *a,int n,int f){
re int l=0,d=1;
for(;d<n;d<<=1)++l;
Rep(i,0,n)
p[i]=(p[i>>1]>>1)^((i&1)<<(l-1));
Rep(i,0,n)
if(i<p[i])
swap(a[i],a[p[i]]);
for(re int i=1;i<n;i<<=1){
cp gn=cp(cos(1.*pi/i),sin(f*1.*pi/i)),w;
for(re int j=0;w=cp(1,0),j<n;j+=(i<<1))
for(re int k=j;k<i+j;++k,w=w*gn){
cp x=a[k],y=w*a[i+k];
a[k]=x+y,a[i+k]=x-y;
}
}
if(f==-1)
Rep(i,0,n)
a[i].r=1.*a[i].r/n;
}
inline void Mul(int *a,int *b,int *c,int len){
re int sqr=sqrt(mod);
re cp temp;
Rep(i,0,len){
a1[i]=a[i]/sqr,b1[i]=a[i]%sqr;
a2[i]=b[i]/sqr,b2[i]=b[i]%sqr;
}
fft(a1,len,1),fft(b1,len,1),fft(a2,len,1),fft(b2,len,1);
Rep(i,0,len){
c1[i]=a1[i]*a2[i];
c2[i]=a1[i]*b2[i]+a2[i]*b1[i];
c3[i]=b1[i]*b2[i];
}
fft(c1,len,-1),fft(c2,len,-1),fft(c3,len,-1);
Rep(i,0,len)
c[i]=((ll)(round(c1[i].r))%mod*sqr%mod*sqr%mod+(ll)(round(c2[i].r))%mod*sqr%mod+(ll)(round(c3[i].r))%mod)%mod;
}
inline void Det(int *a,int *b,int len){
b[len-1]=0;
Rep(i,1,len)
b[i-1]=1ll*a[i]*i%mod;
}
inline void Area(int *a,int *b,int len){
b[0]=0;
Rep(i,1,len)
b[i]=1ll*a[i-1]*inv[i]%mod;
}
inline void Inv(int *a,int *b,int len){
if(len==1){
b[0]=fp(a[0],mod-2);
return;
}
Inv(a,b,len>>1);
Rep(i,0,len)d[i]=a[i],e[i]=b[i];
/*
fft(d,tmp,1),fft(e,tmp,1);
Rep(i,0,tmp)
d[i]=1ll*d[i]*e[i]%mod*e[i]%mod;
fft(d,tmp,-1);
*/
re int tmp=len<<1;
Mul(d,e,d,tmp),Mul(d,e,d,tmp);
Rep(i,0,len)
b[i]=inc(b[i],b[i]),
b[i]=inc(b[i],mod-d[i]);
Rep(i,0,tmp)
d[i]=e[i]=0;
}
inline void Ln(int *a,int *b,int len){
Inv(a,f,len),Det(a,d,len);
re int tmp=len<<1;
/*
fft(f,tmp,1),fft(d,tmp,1);
Rep(i,0,tmp)
f[i]=1ll*f[i]*d[i]%mod;
fft(f,tmp,-1),
*/
Mul(f,d,f,tmp);
Area(f,b,len);
Rep(i,0,tmp)
f[i]=d[i]=0;
}
inline void Exp(int *a,int *b,int len){
if(len==1){
b[0]=1;
return ;
}
Exp(a,b,len>>1),Ln(b,c,len);
Rep(i,0,len)
c[i]=inc(a[i],mod-c[i]),f[i]=b[i];
c[0]=inc(c[0],1);
re int tmp=len<<1;
/*
fft(c,tmp,1),fft(f,tmp,1);
Rep(i,0,tmp)c[i]=1ll*c[i]*f[i]%mod;
fft(c,tmp,-1);
*/
Mul(c,f,c,tmp);
Rep(i,0,len)
b[i]=c[i];
Rep(i,0,tmp)
c[i]=f[i]=0;
}
int main(){
gii(n,m);
inv[1]=1;
rep(i,2,1e5)
inv[i]=mod-1ll*mod/i*inv[mod%i]%mod;
re int d=1,v,ans=0;
rep(i,1,n){
gi(v);
for(re int j=v,k=1;j<=m;j+=v,++k)
a[j]=inc(a[j],inv[k]);
}
while(d<=m)d<<=1;
Exp(a,b,d);
rep(i,1,m)
ans=inc(ans,b[i]);
printf("%d\n",ans);
return 0;
}
code

posted @ 2018-08-13 16:30  Stump  阅读(75)  评论(0编辑  收藏