题目大意：
设$S(n,m)$为第二类斯特林数，$F_i$表示斐波那契数列第$i$项。
给定$n,R,K$,求$sumlimits_{i=1}^{n}(sumlimits_{m=1}^{R}F_i)!i!sumlimits_{l=0}^{i}sumlimits_{j=0}^{sumlimits_{t=1}^{R}F_t}frac{S(k,i-l)}{l!}frac{S(i,sumlimits_{w=1}^{R}F_w-j)}{j!}$的值$mod$ $1000000007$。
其中$n,Rleq10^{18}$，$kleq2*10^5$
题解：
（爽感）
这道题...精神污染。。。当然也包括题解的$O(k)$做法。。。
于是我就自己YY了一个$O(klogk)$的能过的做法.
首先，需要知道三个小结论：
1.$sumlimits_{i=1}^{n}F_i=F_{n+2}-1$，用归纳法可以轻松证明。
2.$m^n=sumlimits_{i=1}^{m}S(n,i)C_m^ii!$
3.$S(n,m)=frac{1}{m!}sumlimits_{i=0}^m(-1)^iC_m^i(m-i)^n$，本质上为结论2经二项式反演后的结果。
证明一下第二个结论：
考虑构造，把$n$个不同的球放入$m$个不同的盒子里，允许有空盒的方案数。
显然，每个球都有$m$种选择方法，答案是$m^n$。
换一种思维方式，可以想成枚举这$n$个球放入了那些盒子里，由于$S(n,m)$表示$n$个不同的球放入$m$个相同的盒子里无空盒的方案数，所以在乘上$C_m^ii!$。
然后就可以开心的化式子了：
设$L=sumlimits_{m=1}^{R}F_i$，再把原式挪一下就变成了：
$sumlimits_{i=1}^{n}sumlimits_{l=0}^{i}i!frac{S(k,i-l)}{l!}sumlimits_{j=0}^LL!frac{S(i,L-j)}{j!}$
然后我们神奇的发现后面的两个式子形式相同，以后一个为例：
$sumlimits_{j=0}^LL!frac{S(i,L-j)}{j!}$，换成枚举$L-j$
$=sumlimits_{j=0}^LL!frac{S(i,j)}{(L-j)!}$
$=sumlimits_{j=0}^LS(i,j)C_L^jj!$，这不就是结论2吗，直接等于$L^i$
因此原式直接变为$sumlimits_{i=0}^nL^ii^k$，但$n$依然很大。
注意当$L=1$时式子变为$sumlimits_{i=0}^ni^k$，需要用拉格朗日插值法求解，详见CF622F.
若$L eq1$，我的做法是把$i^k$变回斯特林数，即$sumlimits_{i=0}^nL^isumlimits_{j=0}^kS(k,j)C_i^jj!$
交换求和号，即$sumlimits_{j=0}^kS(k,j)j!sumlimits_{i=j}^nC_i^jL^i$，考虑如何快速求出$sumlimits_{i=j}^nC_i^jL^i$。
设$g(x)=sumlimits_{i=x}^nC_i^xL^i$
推导一波：
$g(x)=sumlimits_{i=x}^nL^iC_i^x$
$=sumlimits_{i=x}^nL^iC_{i-1}^{x-1}+sumlimits_{i=x}^nL^iC_{i-1}^x$
$=Lsumlimits_{i=x-1}^{n-1}L^{i}C_{i}^{x-1}+Lsumlimits_{i=x}^{n-1}L^{i}C_{i}^x$
$=L(g(x-1)-L^nC_n^{x-1})+L(g(x)-L^nC_n^x)$
在整理一下就是：$g(x)=frac{Lg(x-1)-L^{n+1}C_{n+1}^x}{1-L}$，那么就可以$O(k)$的复杂度求出$g(x)$
到此为止，原式变为$sumlimits_{j=1}^kS(k,j)j!g(j)$，问题就剩下如何求第二类斯特林数$S(k,j)$了。
应用第三个结论：$S(n,m)=frac{1}{m!}sumlimits_{i=1}^m(-1)^iC_m^i(m-i)^n$，将其化为卷积形式：
$S(n,m)=sumlimits_{i=1}^mfrac{(-1)^i}{i!}frac{(m-i)^n}{(m-i)!}$$NTT$即可。
又由于模数不是$NTT$模数，因此需要$MTT$来实现。
#include<cmath> #include<cstdio> #include<algorithm> #define N 530010 #define Re register #define ld long double #define mod 1000000007 #define pi (ld)acos(-1) #pragma GCC optimize(3) #pragma GCC optimize("inline") #pragma GCC optimize("Ofast") using namespace std; typedef long long ll; struct node{ll a[2][2];}tt,an; int k,R[N]; ll n,r,ans,y[N],S[N],fac[N],inv[N]; node mul(const node&x,const node&y) { node ret; for(Re int i = 0;i<=1;i++) { for(Re int j = 0;j<=1;j++) { ret.a[i][j] = 0; for(Re int k = 0;k<=1;k++) ret.a[i][j] = (ret.a[i][j]+x.a[i][k]*y.a[k][j]%mod)%mod; } }return ret; } inline void mpow(ll y) { an.a[0][0] = an.a[1][1] = 1; while(y) { if(y&1)an = mul(an,tt); tt = mul(tt,tt),y>>=1; } } ll ksm(ll x,ll y) { ll fin = 1; x%=mod,y%=mod-1; while(y) { if(y&1)fin = fin*x%mod; x = x*x%mod,y>>=1; }return fin; } struct fs { ld x,y; friend fs operator+(const fs&a,const fs&b){return (fs){a.x+b.x,a.y+b.y};} friend fs operator-(const fs&a,const fs&b){return (fs){a.x-b.x,a.y-b.y};} friend fs operator*(const fs&a,const fs&b){return (fs){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};} }a[N],b[N],c[N],d[N],e[N],f[N],g[N],h[N]; inline void FFT(fs*A,int lim,int fl) { for(Re int i = 0;i<lim;i++)if(i<R[i])swap(A[i],A[R[i]]); for(Re int i = 2;i<=lim;i<<=1) { fs wn = (fs){cos(2*pi/i),fl*sin(2*pi/i)}; for(Re int j = 0;j<lim;j+=i) { fs w = (fs){1,0}; for(Re int k = 0;k<i>>1;k++) { fs x = A[j+k],y = w*A[j+k+(i>>1)]; A[j+k] = x+y,A[j+k+(i>>1)] = x-y; w = w*wn; } } }if(fl==-1)for(Re int i = 0;i<lim;i++)A[i].x/=lim; } inline void p2() { for(Re int i = 1;i<=k+2;i++)y[i] = (y[i-1]+ksm(i,k))%mod; for(Re int i = 1;i<=k+2;i++) { int t = 1; for(Re int j = 1;j<=k+2;j++) { if(i==j)continue; t = t*(n%mod-j)%mod*ksm(i-j,mod-2)%mod; }ans = (ans+y[i]*t%mod)%mod; }printf("%lld",(ans+mod)%mod),exit(0); } int main() { freopen("c.in","r",stdin); freopen("c.out","w",stdout); scanf("%lld%lld%d",&n,&r,&k); tt.a[0][0] = tt.a[0][1] = tt.a[1][0] = 1; mpow(r+1); r = an.a[0][0]-1; if(r==1)p2(); fac[0] = fac[1] = inv[0] = inv[1] = 1; for(Re int i = 2;i<=k;i++)fac[i] = fac[i-1]*i%mod,inv[i] = (mod-mod/i)*inv[mod%i]%mod; for(Re int i = 2;i<=k;i++)inv[i] = inv[i]*inv[i-1]%mod; for(Re int i = 0;i<=k;i++) { ll A = (inv[i]*(i&1?-1:1)+mod)%mod; ll B = ksm(i,k)*inv[i]%mod; a[i].x = (ld)(A>>15); b[i].x = (ld)(A&0x7fff); c[i].x = (ld)(B>>15); d[i].x = (ld)(B&0x7fff); } int lim = 1,bit = 0; while(lim<=k<<1)lim<<=1,bit++; for(Re int i = 0;i<lim;i++)R[i] = (R[i>>1]>>1)|((i&1)<<bit-1); FFT(a,lim,1),FFT(b,lim,1),FFT(c,lim,1),FFT(d,lim,1); for(Re int i = 0;i<lim;i++) { e[i] = a[i]*c[i]; f[i] = a[i]*d[i]+b[i]*c[i]; g[i] = b[i]*d[i]; } FFT(e,lim,-1),FFT(f,lim,-1),FFT(g,lim,-1); for(Re int i = 0;i<=k;i++) { ll E = (ll)round(e[i].x)%mod,F = (ll)round(f[i].x)%mod,G = (ll)round(g[i].x)%mod; S[i] = (((E<<30)%mod+(F<<15)%mod+G)%mod+mod)%mod; } ll lst = (ksm(r,n+1)-1)*ksm(r-1,mod-2)%mod,now = 1,t = ksm(r,n+1),tt = ksm(1-r,mod-2); for(Re int i = 1;i<=k;i++) { now = now*(n%mod+2-i)%mod; lst = (lst*r%mod-now*inv[i]%mod*t%mod)%mod*tt%mod; ans = (ans+S[i]*fac[i]%mod*lst%mod)%mod; }printf("%lld",(ans+mod)%mod); return 0; }