51nod1229 序列求和 V2

Posted 2021-02-17 judge

这题...毒瘤吧，可能要写两份代码...

传送门

noteskey

我们考虑这里的复杂度肯定是与 k 相关的，而且平方也是没问题的，那么我们先看看 S(k) 能怎么得到：

[egin{aligned}S(k)=&sum_{i=1}^n i^k r^i\ r·S(k)=&sum_{i=2}^{n+1} (i-1)^k r^i \ (r-1)S(k)=& r^{n+1}n^k-r+ sum_{j=2}^{k-1} r^i((j-1)^k-j^k)\\ & ext{然后二项式展开}~ (j-1)^k ~ ext{可得} \\ (r-1)·S(k)=& r^{n+1} n^k-r+sum_{j=0}^{k-1} (-1)^{k-j} egin{pmatrix} k\ j end{pmatrix} Big(sum_{i=2}^{n}i^jr^i Big) \\ & ext{不难发现后面那个玩意儿其实是} ~S(j)-r \\ (r-1)·S(k)=& r^{n+1} n^k-r+sum_{j=0}^{k-1} (-1)^{k-j} egin{pmatrix} k \ j end{pmatrix} ig( S(j)-rig) \ S(k)=& {r^{n+1} n^k-r+sum_{j=0}^{k-1} (-1)^{k-j} egin{pmatrix} k \ j end{pmatrix} ig( S(j)-rig) over r-1} end{aligned} ]

这样的话我们就可以递推出解了...复杂度 (n^2) 很 ok

慢着！这样真的不锅嘛？ r=1 怎么办...

对哟，所以我们还要打一份 (O(n)) 的拉格朗日插值代码...来特判 r=1 的情况，因为 r=1 的时候问题就变成了求幂和函数 (S(k)) 的值

当然你也可以选择打个斯特林数的 (n^2) 代码哟，不过想用伯努利的同学我真的不想说什么...您看看这题的模数丫，莫非您还要来个 MTT 然后复杂度 还比不过 拉格朗日的说？

还有 r=0 的情况？ 这个时候答案必然为 0 算个鸡儿...

code

巨丑无比...

//by Judge
#pragma GCC optimize("Ofast")
#include<cstdio>
#include<cstring>
#include<iostream>
#define Rg register
#define fp(i,a,b) for(Rg int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(Rg int i=(a),I=(b)-1;i>I;--i)
#define ll long long
using namespace std;
const int mod=1e9+7;
const int M=2003;
typedef int arr[M];
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
char buf[1<<21],*p1=buf,*p2=buf;
inline int MO(Rg ll x){return x-x/mod*mod;}
inline int mul(Rg int x,Rg int y){return MO(1ll*x*y);}
inline int dec(Rg int x,Rg int y){return x<y?x-y+mod:x-y;}
inline int inc(Rg int x,Rg int y){return x+y>=mod?x+y-mod:x+y;}
inline ll read(){ ll x=0; char c=getchar(); while(!isdigit(c)) c=getchar();
    for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x;
} char sr[1<<21],z[20];int CCF=-1,Z;
inline void Ot(){fwrite(sr,1,CCF+1,stdout),CCF=-1;}
inline void print(int x,char chr='
'){ if(CCF>1<<20)Ot();
    while(z[++Z]=x%10+48,x/=10); while(sr[++CCF]=z[Z],--Z);sr[++CCF]=chr;
} ll n,r; int w,k,ans; arr s,fac,ifac,pre,suf;
inline int qpow(Rg int x,Rg int p=mod-2,Rg int s=1){
    for(;p;p>>=1,x=mul(x,x)) if(p&1) s=mul(s,x); return s;
}
inline int C(int n,int m){return mul(fac[n],mul(ifac[m],ifac[n-m]));}
inline int calc(){ s[0]=0,n%=mod; //拉格朗日插值
    fp(i,1,k+2) s[i]=inc(s[i-1],qpow(i,k));
    pre[0]=suf[k+3]=w=1,ans=0;
    fp(i,1,k+2) pre[i]=mul(pre[i-1],dec(n,i));
    fd(i,k+2,1) suf[i]=mul(suf[i+1],dec(n,i));
    fp(i,2,k+2) w=mul(w,qpow(dec(1,i)));
    fp(i,1,k+2) ans=inc(ans,mul(mul(s[i],w),mul(pre[i-1],suf[i+1]))),
        w=mul(mul(w,dec(i,inc(k,2))),mul(ifac[i],fac[i-1])); return ans;
}
inline int solv(){ if(!r) return 0; // n^2 递推
    Rg int m=(n+1>=mod-1)?(n+1)%(mod-1)+mod-1:n+1;
    s[0]=mul(dec(qpow(r,m),r),qpow(r-1));
    fp(i,1,k){ s[i]=dec(mul(qpow(r,m),qpow(n%mod,i)),r);
        fp(j,0,i-1) if((i^j)&1)
            s[i]=dec(s[i],mul(C(i,j),dec(s[j],r)));
            else s[i]=inc(s[i],mul(C(i,j),dec(s[j],r)));
        s[i]=mul(s[i],qpow(r-1));
    } return s[k];
}
inline void prep(int n){ fac[0]=fac[1]=ifac[0]=ifac[1]=1;
    fp(i,2,n) ifac[i]=mul(mod-mod/i,ifac[mod%i]);
    fp(i,2,n) fac[i]=mul(fac[i-1],i),ifac[i]=mul(ifac[i],ifac[i-1]);
}
int main(){ prep(2000); int T=read();
    for(;T;--T) n=read(),k=read(),r=read()%mod,
        print(r^1?solv():calc());
    return Ot(),0;
}