代码

#include<bits/stdc++.h>
#define fo(i,a,b)for(int i=a,_e=b;i<=_e;++i)
#define fd(i,a,b)for(int i=b,_e=a;i>=_e;--i)
#define ff(i,a,b)for(int i=a,_e=b;i<_e;++i)
#define cl(c,Q) memset(c,0,Q*4)
#define ll long long
#define db double
using namespace std;
const int N=(1<<20)+5,mo=1e9+7;
const db pi=acos(-1);
struct Z{
	db x,y;
	Z(db _x=0,db _y=0){x=_x;y=_y;}
	Z operator+(const Z &b){return Z(x+b.x,y+b.y);}
	Z operator-(const Z &b){return Z(x-b.x,y-b.y);}
	Z operator*(const Z &b){return Z(x*b.x-y*b.y,x*b.y+y*b.x);}
}w[N],a[N];
Z conj(const Z &b){return Z(b.x,-b.y);}
int n,k,iv[N];
int h[N],Q;
int f[N],f2[N];
void dft(Z *a,int si){
	ff(i,1,Q)if(h[i]>i)swap(a[i],a[h[i]]);Z A;
	for(int i=1;i<Q;i<<=1)for(int j=0;j<Q;j+=i*2)
		ff(k,0,i)A=w[i+k]*a[i+j+k],a[i+j+k]=a[j+k]-A,a[j+k]=a[j+k]+A;
	if(si==-1){
		reverse(a+1,a+Q);db y=1./Q;
		ff(i,0,Q)a[i].x*=y,a[i].y*=y;
	}
}
void pre(int n){
	for(Q=1;Q<=n;Q<<=1);
	fo(i,1,Q)h[i]=(h[i>>1]>>1)|(i&1?Q>>1:0);
}
void fft(int *a,Z *a0,Z *a1){
	static Z aa[N];
	ff(i,0,Q)aa[i]=Z(a[i]&32767,a[i]>>15);
	dft(aa,1);
	ff(i,0,Q){
		int j=(Q-1)&(Q-i);
		a0[i]=(aa[i]+conj(aa[j]))*Z(0.5,0);
		a1[i]=(aa[i]-a0[i])*Z(0,-1);
	}
}
void dot(int *c,Z *a0,Z *a1,Z *b0,Z *b1){
	static Z a,b;
	ff(i,0,Q){
		a=a0[i]*b0[i]+a1[i]*b1[i]*Z(0,1);
		b=a0[i]*b1[i]+a1[i]*b0[i];
		a0[i]=a;a1[i]=b;
	}
	dft(a0,-1);dft(a1,-1);
	ff(i,0,Q)
		c[i]=((ll)(a0[i].x+0.5)+((ll)(a1[i].x+0.5)%mo<<15)+((ll)(a0[i].y+0.5)%mo<<30))%mo;
}
void mul(int *c,int *a,int *b){
	static Z a0[N],a1[N],b0[N],b1[N];
	fft(a,a0,a1);
	fft(b,b0,b1);
	dot(c,a0,a1,b0,b1);
}
void exp_inv(int *b,int *a,int E){
	static int c[N],g[N],p[N];
	if(E==2){b[0]=p[0]=1;return;}
	static Z c0[N],c1[N],g0[N],g1[N];
	ff(i,0,E/4)b[i]=p[i];
	for(int O=E/2;O<=E;O<<=1){
		pre(O-1);
		ff(i,0,Q)c[i]=a[i];
		ff(i,0,Q/2)g[i]=b[i];
		ff(i,Q/2,Q)g[i]=0;
		fft(c,c0,c1);
		fft(g,g0,g1);
		// warning g0,g1
		dot(c,c0,c1,g0,g1);
		ff(i,0,Q/2)c[i]=0;
		fft(c,c0,c1);
		dot(c,c0,c1,g0,g1);
		ff(i,Q/2,Q)b[i]=mo-c[i];
		if(O*2==E)ff(i,0,Q)p[i]=b[i];
	}
}
void exp(int *b,int *a){
	static int c[N],g[N];
	b[0]=1;int E=Q;
	for(int O=2;O<=E;O<<=1){
		exp_inv(g,b,O);
		pre(O-1);
		cl(c,Q);
		ff(i,0,Q/2)c[i]=(ll)b[i+1]*(i+1)%mo;
		mul(c,c,g);
		fd(i,0,Q-2)c[i+1]=(ll)c[i]*iv[i+1]%mo;
		ff(i,Q/2,Q)c[i]=((ll)a[i]-c[i]+mo)%mo;
		cl(c,Q/2);cl(g,Q);
		ff(i,0,Q/2)g[i]=b[i];
		mul(c,c,g);
		ff(i,Q/2,Q)b[i]=c[i];
	}
}
int main(){
	//scanf("%d%d",&n,&k);
	k=n=1e6;
	for(Q=1;Q<=k;Q<<=1);
	for(int i=1;i<Q;i<<=1)
		fo(j,0,i-1)w[i+j]=Z(cos(pi*j/i),sin(pi*j/i));
	
	iv[1]=1;
	fo(i,2,Q)iv[i]=mo-(ll)iv[mo%i]*(mo/i)%mo;
	fo(i,1,k)f[i]=(ll)iv[i]*(n-1)%mo;
	fo(i,2,n)
		fo(j,1,k/i)
			f[i*j]=(f[i*j]+mo-iv[j])%mo;
	
	exp(f2,f);
	printf("%d\n",f2[k]);
}

评测结果