【51Nod1258】序列求和V4（FFT）

题面

51Nod

多组数据，求：

\[Ans=\sum_{i=1}^ni^k,n\le 10^{18},k\le50000
\]

题解

预处理伯努利数，时间复杂度\(O(nlogn)\)

然后利用伯努利数求和即可。

\[\sum_{i=1}^n i^k=\frac{1}{k+1}\sum_{i=0}^kB_iC_{k+1}^i(n+1)^{k+1-i}
\]

预处理需要多项式求逆，因为模数不太好，所以需要\(MTT\)

#include<iostream>

#include<cstdio>

#include<cstdlib>

#include<cstring>

#include<cmath>

#include<algorithm>

#include<set>

#include<map>

#include<vector>

#include<queue>

using namespace std;

#define ll long long

#define RG register

#define MOD 1000000007

#define MAX 150000

const int NN=50000;

const int M=sqrt(MOD);

const double Pi=acos(-1);

inline ll read()

{

    RG ll x=0,t=1;RG char ch=getchar();

    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();

    if(ch=='-')t=-1,ch=getchar();

    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();

    return x*t;

}

int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}

struct Complex{double a,b;}W[MAX],A1[MAX],A2[MAX],B1[MAX],B2[MAX],X[MAX],Y[MAX],Z[MAX];

Complex operator+(Complex a,Complex b){return (Complex){a.a+b.a,a.b+b.b};}

Complex operator-(Complex a,Complex b){return (Complex){a.a-b.a,a.b-b.b};}

Complex operator*(Complex a,Complex b){return (Complex){a.a*b.a-a.b*b.b,a.b*b.a+a.a*b.b};}

int r[MAX],N,l;

void FFT(Complex *P,int N,int opt)

{

	for(int i=0;i<N;++i)if(i<r[i])swap(P[i],P[r[i]]);

	for(int i=1;i<N;i<<=1)

		for(int p=i<<1,j=0;j<N;j+=p)

			for(int k=0;k<i;++k)

			{

				Complex w=(Complex){W[N/i*k].a,W[N/i*k].b*opt};

				Complex X=P[j+k],Y=P[i+j+k]*w;

				P[j+k]=X+Y;P[i+j+k]=X-Y;

			}

	if(opt==-1)for(int i=0;i<N;++i)P[i].a/=N;

}

void MTT(int *a,int *b,int len,int *c)

{

	memset(A1,0,sizeof(A1));memset(B1,0,sizeof(B1));

	memset(A2,0,sizeof(A2));memset(B2,0,sizeof(B2));

	for(N=1,l=0;N<=len;N<<=1)++l;

	for(int i=0;i<N;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));

	for(int i=1;i<N;i<<=1)

		for(int k=0;k<i;++k)W[N/i*k]=(Complex){cos(k*Pi/i),sin(k*Pi/i)};

	for(int i=0;i<len;++i)a[i]%=MOD,b[i]%=MOD;

	for(int i=0;i<len;++i)A1[i].a=a[i]/M,A2[i].a=a[i]%M;

	for(int i=0;i<len;++i)B1[i].a=b[i]/M,B2[i].a=b[i]%M;

	memset(X,0,sizeof(X));memset(Y,0,sizeof(Y));memset(Z,0,sizeof(Z));

	FFT(A1,N,1);FFT(A2,N,1);FFT(B1,N,1);FFT(B2,N,1);

	for(int i=0;i<N;++i)

	{

		X[i]=A1[i]*B1[i];

		Y[i]=A1[i]*B2[i]+A2[i]*B1[i];

		Z[i]=A2[i]*B2[i];

	}

	FFT(X,N,-1);FFT(Y,N,-1);FFT(Z,N,-1);

	for(int i=0;i<len;++i)

	{

		int ans=0;

		ans=(ll)(X[i].a+0.5)%MOD*M%MOD*M%MOD;

		ans=(ans+(ll)(Y[i].a+0.5)%MOD*M)%MOD;

		ans=(ans+(ll)(Z[i].a+0.5))%MOD;

		c[i]=ans;

	}

}

int c[MAX],d[MAX];

void Inv(int *a,int *b,int len)

{

	if(len==1){b[0]=fpow(a[0],MOD-2);return;}

	Inv(a,b,len>>1);

	MTT(a,b,len,c);MTT(b,c,len,d);

	for(int i=0;i<len;++i)b[i]=(b[i]+b[i])%MOD;

	for(int i=0;i<len;++i)b[i]=(b[i]+MOD-d[i])%MOD;

}

int a[MAX];

int n,B[MAX],jc[MAX],jv[MAX],inv[MAX];

int C(int n,int m){return 1ll*jc[n]*jv[m]%MOD*jv[n-m]%MOD;}

int main()

{

	B[0]=jc[0]=jv[0]=inv[0]=inv[1]=1;

	for(int i=2;i<=NN+NN;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;

	for(int i=1;i<=NN+NN;++i)jc[i]=1ll*jc[i-1]*i%MOD;

	for(int i=1;i<=NN+NN;++i)jv[i]=1ll*jv[i-1]*inv[i]%MOD;

	for(int i=0;i<=NN;++i)a[i]=jv[i+1];

	Inv(a,B,1<<16);

	for(int i=0;i<=NN;++i)B[i]=1ll*B[i]*jc[i]%MOD;

	int T=read();

	while(T--)

	{

		int n=read()%MOD,k=read(),ans=0,nw=n+1;

		for(int i=k;~i;--i,nw=1ll*nw*(n+1)%MOD)ans=(ans+1ll*C(k+1,i)*B[i]%MOD*nw)%MOD;

		ans=1ll*ans*inv[k+1]%MOD;

		printf("%d\n",ans);

	}

	return 0;

}

巴特西

【51Nod1258】序列求和V4（FFT）

【51Nod1258】序列求和V4（FFT）

题面

题解

最新文章

热门文章