首页 > 代码库 > 51nod1258 序列求和V4
51nod1258 序列求和V4
T(n) = n^k，S(n) = T(1) + T(2) + ...... T(n)。给出n和k，求S(n)。
例如k = 2，n = 5，S(n) = 1^2 + 2^2 + 3^2 + 4^2 + 5^2 = 55。
由于结果很大，输出S(n) Mod 1000000007的结果即可。
Input
第1行：一个数T，表示后面用作输入测试的数的数量。（1 <= T <= 500)第2 - T + 1行：每行2个数，N, K中间用空格分割。(1 <= N <= 10^18, 1 <= K <= 50000)Output
共T行，对应S(n) Mod 1000000007的结果。设伯努利数第i项为B[i]，则有结论
由于伯努利数的指数生成函数为
，化简得
，计算其在模x50001意义下，系数模1e9+7的结果即得到B[0..50000] mod (1e9+7)
由于模数较大，计算多项式逆元可以用倍增+（三模数ntt+CRT合并（每次倍增需15次ntt）或分段fft（一般写法每次倍增需12次dft）），复杂度为O(50000log50000)，常数较大
#include<cstdio>#include<cmath>#include<algorithm>typedef long double ld;typedef long long i64;const ld pi=3.1415926535897932384626433832795l;const int P=1000000007,M=31623;struct cplx{ ld a,b; cplx(ld _a=0.l,ld _b=0.l):a(_a),b(_b){} cplx operator+(cplx x){return cplx(a+x.a,b+x.b);} cplx operator-(cplx x){return cplx(a-x.a,b-x.b);} cplx operator*(cplx x){return cplx(a*x.a-b*x.b,a*x.b+b*x.a);} cplx operator~(){return cplx(a,-b);}}w[2][131072],A[131072],B[131072],C[131072],D[131072];int N,X,r[131072];void dft(cplx*a,int t){ for(int i=0;i<N;++i)if(i>r[i])std::swap(a[i],a[r[i]]); for(int i=1;i<N;i<<=1){ for(int j=0,v=65536/i;j<N;j+=i<<1){ cplx*b=a+j,*c=b+i; for(int k=0,p=0;k<i;++k,p+=v){ cplx x=b[k],y=c[k]*w[t][p]; b[k]=x+y; c[k]=x-y; } } } if(t)for(int i=0;i<N;++i)a[i].a/=N;}int pow(int a,int n){ int v=1; for(;n;n>>=1){ if(n&1)v=i64(v)*a%P; a=i64(a)*a%P; } return v;}int fac[50007]={1},fiv[50007]={1},v1[50007],v2[50007],v3[50007],b[50007];void calc(int n){ if(n==1){ v2[0]=1; return; } int n1=n+1>>1; calc(n1); for(N=2,X=0;N<n*2;N<<=1,++X); for(int i=1;i<N;++i)r[i]=r[i>>1]>>1|(i&1)<<X; for(int i=0;i<n1;++i)A[i]=cplx(v2[i]/M),B[i]=cplx(v2[i]%M); for(int i=n1;i<N;++i)A[i]=B[i]=cplx(); dft(A,0);dft(B,0); for(int i=0;i<N;++i){ C[i]=B[i]*B[i]; B[i]=A[i]*B[i]; A[i]=A[i]*A[i]; } dft(A,1);dft(B,1);dft(C,1); for(int i=0;i<n;++i)v3[i]=(i64(A[i].a+0.4999l)%P*(M*M)%P+i64(B[i].a+0.4999l)%P*2*M%P+i64(C[i].a+0.4999l))%P; for(int i=0;i<n;++i)A[i]=cplx(v3[i]/M),B[i]=cplx(v3[i]%M),C[i]=cplx(v1[i]/M),D[i]=cplx(v1[i]%M); for(int i=n;i<N;++i)A[i]=B[i]=C[i]=D[i]=cplx(); dft(A,0);dft(B,0);dft(C,0);dft(D,0); for(int i=0;i<N;++i){ cplx x=A[i]*D[i]+B[i]*C[i]; A[i]=A[i]*C[i]; C[i]=B[i]*D[i]; B[i]=x; } dft(A,1);dft(B,1);dft(C,1); for(int i=0;i<n;++i){ int x=(v2[i]*2ll+P-(i64(A[i].a+0.4999l)%P*(M*M)%P+i64(B[i].a+0.4999l)%P*M%P+i64(C[i].a+0.4999l)%P))%P; v2[i]=x<0?x+P:x; }}int _C(int n,int m){ return fac[n]*1ll*fiv[m]%P*fiv[n-m]%P;}int main(){ for(int i=0;i<131072;++i){ w[0][i]=cplx(cos(pi*i/65536.l),sin(pi*i/65536.l)); w[1][i]=~w[0][i]; } for(N=2,X=0;N<111;N<<=1,++X); for(int i=1;i<N;++i)r[i]=r[i>>1]>>1|(i&1)<<X; for(int i=1;i<=50001;++i){ fac[i]=1ll*fac[i-1]*i%P; fiv[i]=pow(fac[i],P-2); } for(int i=0;i<=50000;++i)v1[i]=fiv[i+1]; calc(50001); for(int i=0;i<=50000;++i)b[i]=v2[i]*1ll*fac[i]%P; i64 n0,n;int T,k; for(scanf("%d",&T);T;--T){ scanf("%lld%d",&n0,&k); int ans=0; n=(n0+1)%P; for(int i=1,pw=n;i<=k+1;++i,pw=1ll*pw*n%P){ ans=(ans+1ll*_C(k+1,i)*b[k+1-i]%P*pw)%P; } ans=1ll*ans*pow(k+1,P-2)%P; if(ans<0)ans+=P; printf("%d\n",ans); } return 0;}
51nod1258 序列求和V4
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。