[51nod 1258] [伯努利數] [多項式求逆] [任意模數NTT] 序列求和 V4
阿新 • • 發佈:2019-01-31
上次做一套模擬賽的時候,其中需要求自然數k次冪和,然後我只會n^2的…我記得n^2有20分,nlogn求可以爆到90分……
——鏼鏼鏼2015年國家集訓隊論文
大概就這樣多項式求個逆,求出生成函式就可以了
這樣是
至於為什麼要大於
不過直接上CRT要爆longlong,可以用一個小技巧處理
這樣每次NTT分成三個素數的NTT,再合併,就可以了
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <string>
#include <cstring>
using namespace std;
typedef long long ll;
const int N=300010,P=1e9+7;
const int p1=998244353,p2=1004535809,p3=469762049 ;
inline char nc(){
static char buf[100000],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}
inline void rea(int &x){
char c=nc(); x=0;
for(;c>'9'||c<'0';c=nc());for(;c>='0'&&c<='9';x=x*10+c-'0',c=nc());
}
inline void rea(ll &x ){
char c=nc(); x=0;
for(;c>'9'||c<'0';c=nc());for(;c>='0'&&c<='9';x=x*10+c-'0',c=nc());
}
inline ll Pow(ll x,ll y,ll p){
x%=p; ll ret=1;
for(;y;y>>=1,x=x*x%p) if(y&1) ret=ret*x%p;
return ret;
}
const ll M_=1LL*p1*p2,inv1=Pow(p1,p2-2,p2),inv2=Pow(p2,p1-2,p1),inv3=Pow(M_%p3,p3-2,p3);
int rev[N],fac[N],inv[N],a[N],b[N];
struct NTT{
int P,num;
int w[2][N];
void Pre(int p,int n){
P=p; num=n;
int ng=Pow(3,(P-1)/num,P),ig=Pow(ng,P-2,P);
w[1][0]=w[0][0]=1;
for(int i=1;i<num;i++)
w[1][i]=1LL*w[1][i-1]*ng%P,w[0][i]=1LL*w[0][i-1]*ig%P;
}
void FFT(int *a,int n,int r){
for(int i=1;i<n;i++) if(rev[i]>i) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1)
for(int j=0;j<n;j+=(i<<1))
for(int k=0;k<i;k++){
int x=a[j+k],y=1LL*a[j+k+i]*w[r][num/(i<<1)*k]%P;
a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P;
}
if(!r) for(int i=0,inv=Pow(n,P-2,P);i<n;i++) a[i]=1LL*a[i]*inv%P;
}
}ntt[3];
inline void Pre(int n){
fac[0]=1; for(int i=1;i<=n;i++) fac[i]=1LL*fac[i-1]*i%P;
inv[1]=1; for(int i=2;i<=n;i++) inv[i]=1LL*(P-P/i)*inv[P%i]%P;
inv[0]=1; for(int i=1;i<=n;i++) inv[i]=1LL*inv[i]*inv[i-1]%P;
}
inline ll mul(ll a,ll b,ll p){
a%=p; b%=p;
return ((a*b-(ll)((ll)((long double)a/p*b+1e-3)*p))%p+p)%p;
}
inline int CRT(int a1,int a2,int a3){
ll A=(mul(1LL*a1*p2%M_,inv2,M_)+mul(1LL*a2*p1%M_,inv1,M_))%M_;
ll K=(1LL*a3+p3-A%p3)*inv3%p3;
return (K*(M_%P)+A)%P;
}
int _b[3][N],__b[N];
void Inv(int *a,int *b,int n){
static int tmp[N];
if(n==1) return void(b[0]=CRT(Pow(a[0],p1-2,p1),Pow(a[0],p2-2,p2),Pow(a[0],p3-2,p3)));
Inv(a,b,n>>1);
int L=0; while(!(n>>L&1)) L++;
for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
for(int k=0;k<3;k++)
for(int i=0;i<n;i++) _b[k][i]=b[i],_b[k][i+n]=0;
for(int k=0;k<3;k++){
for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0;
ntt[k].FFT(tmp,n<<1,1); ntt[k].FFT(_b[k],n<<1,1);
for(int i=0;i<(n<<1);i++) _b[k][i]=1LL*_b[k][i]*tmp[i]%ntt[k].P;
ntt[k].FFT(_b[k],n<<1,0);
}
for(int i=0;i<n;i++) _b[0][i]=_b[1][i]=_b[2][i]=(P-CRT(_b[0][i],_b[1][i],_b[2][i]))%P;
_b[0][0]=_b[1][0]=_b[2][0]=(_b[0][0]+2)%P;
for(int k=0;k<3;k++){
for(int i=0;i<n;i++) tmp[i]=b[i],tmp[i+n]=0;
ntt[k].FFT(tmp,n<<1,1); ntt[k].FFT(_b[k],n<<1,1);
for(int i=0;i<(n<<1);i++) _b[k][i]=1LL*_b[k][i]*tmp[i]%ntt[k].P;
ntt[k].FFT(_b[k],n<<1,0);
}
for(int i=0;i<n;i++) b[i]=CRT(_b[0][i],_b[1][i],_b[2][i]);
}
int t,k;
ll n;
inline int C(int x,int y){
return 1LL*fac[x]*inv[y]%P*inv[x-y]%P;
}
int main(){
int M=1; for(;M<=50000;M<<=1);
ntt[0].Pre(p1,M<<1);
ntt[1].Pre(p2,M<<1);
ntt[2].Pre(p3,M<<1);
Pre(M);
for(int i=0;i<M;i++) a[i]=inv[i+1];
Inv(a,b,M);
for(int i=0;i<M;i++) b[i]=1LL*b[i]*fac[i]%P;
b[1]=P-b[1];
rea(t);
while(t--){
rea(n); rea(k); n%=P;
int ans=0;
for(int i=k,prod=n;~i;prod=1LL*prod*n%P,i--)
ans=(1LL*prod*b[i]%P*C(k+1,i)%P+ans)%P;
ans=1LL*ans*Pow(k+1,P-2,P)%P;
printf("%d\n",ans);
}
return 0;
}