以前SB了写过一篇比较复杂的做法,现在才发现这个东西非常的Naive。
设前n个选完之后剩下 ,我们要求 。
直接把这个东西用第一类斯特林数搞一搞,我们要求
考虑维护 ,发现前n个每加入一个数,就变为 ,拆开是一个多项式卷积,我们用多项式快速幂即可在 解决(后者先求对数再exp回去即可)。
不过瓶颈在于这个 从 取到 ,还要自然数幂求和,时间复杂度
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1e3+50, mod=1e9+7;
inline int add(int x,int y) {return (x+y>=mod) ? (x+y-mod) : (x+y);}
inline int dec(int x,int y) {return (x-y<0) ? (x-y+mod) : (x-y);}
inline int mul(int x,int y) {return (LL)x*y%mod;}
inline int power(int a,int b,int rs=1) {for(;b;b>>=1,a=mul(a,a)) if(b&1) rs=mul(rs,a); return rs;}
int n,m,t,k; LL s;
int S[N][N],fac[N],ifac[N],inv[N];
struct lgip {
int pre[N],suf[N],yc[N],deg;
inline void init() {
for(int i=1;i<=k+2;i++) yc[i]=add(yc[i-1],yc[i]);
}
inline int getval(int nn) {
int rs=0, o=1;
for(int i=2;i<=k+2;i++) o=mul(o,inv[i-1]);
pre[0]=1; for(int i=1;i<=k+2;i++) pre[i]=mul(pre[i-1],dec(i,nn));
suf[k+3]=1; for(int i=k+2;i;i--) suf[i]=mul(suf[i+1],dec(i,nn));
for(int i=1;i<=k+2;i++) {
rs=add(rs,mul(mul(pre[i-1],suf[i+1]),mul(o,yc[i])));
o=mul(o,k+2-i); o=mul(o,mod-inv[i]);
} return rs;
}
} f[N];
struct poly {
int a[N];
poly() {memset(a,0,sizeof(a));}
friend inline poly operator *(const poly &a,const poly &b) {
poly c;
for(int i=0;i<=k;i++)
for(int j=0;j<=i;j++)
c.a[i]=add(c.a[i],mul(a.a[j],b.a[i-j]));
return c;
}
} F,G,H;
int main() {
cin>>s>>t>>n>>m; k=m-n; s=s%mod;
fac[0]=1; for(int i=1;i<=k+2;i++) fac[i]=mul(fac[i-1],i);
ifac[k+2]=power(fac[k+2],mod-2); for(int i=k+1;~i;i--) ifac[i]=mul(ifac[i+1],i+1);
inv[0]=1; for(int i=1;i<=k+2;i++) inv[i]=mul(fac[i-1],ifac[i]);
for(int i=1;i<=k+2;i++)
for(int pw=1,j=0;j<=k;j++) {
f[j].yc[i]=pw;
pw=mul(pw,i);
}
for(int i=0;i<=k;i++) f[i].init();
for(int i=0;i<=k;i++) {
G.a[i]=mul(ifac[i],f[i].getval(t));
if(i&1) G.a[i]=mod-G.a[i];
} H.a[0]=1;
for(;n;n>>=1,G=G*G)
if(n&1) H=G*H;
for(int i=0;i<=k;i++) F.a[i]=power(s,i);
for(int i=0;i<=k;i++) F.a[i]=mul(F.a[i],ifac[i]);
F=F*H; int rs=0;
for(int i=0;i<=k;i++) F.a[i]=mul(F.a[i],fac[i]);
S[0][0]=1;
for(int i=1;i<=k;i++)
for(int j=1;j<=k;j++)
S[i][j]=add(S[i-1][j-1],mul(S[i-1][j],i-1));
for(int i=0;i<=k;i++) {
if((k-i)&1) rs=dec(rs,mul(S[k][i],F.a[i]));
else rs=add(rs,mul(S[k][i],F.a[i]));
} cout<<mul(rs,ifac[k])<<'\n';
}