2019.10.02模拟赛T3

题目大意：

　　设$S(n,m)$为第二类斯特林数，$F_i$表示斐波那契数列第$i$项。

　　给定$n,R,K$,求$\sum\limits_{i=1}^{n}(\sum\limits_{m=1}^{R}F_i)!i!\sum\limits_{l=0}^{i}\sum\limits_{j=0}^{\sum\limits_{t=1}^{R}F_t}\frac{S(k,i-l)}{l!}\frac{S(i,\sum\limits_{w=1}^{R}F_w-j)}{j!}$的值$mod$ $1000000007$。

　　其中$n,R\leq10^{18}$，$k\leq2*10^5$

题解：

　　（爽感）

　　这道题...精神污染。。。当然也包括题解的$O(k)$做法。。。

　　于是我就自己YY了一个$O(klogk)$的能过的做法.

　　首先，需要知道三个小结论：

　　　　1.$\sum\limits_{i=1}^{n}F_i=F_{n+2}-1$，用归纳法可以轻松证明。

　　　　2.$m^n=\sum\limits_{i=1}^{m}S(n,i)C_m^ii!$

　　　　3.$S(n,m)=\frac{1}{m!}\sum\limits_{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=\sum\limits_{m=1}^{R}F_i$，再把原式挪一下就变成了：

　　　　$\sum\limits_{i=1}^{n}\sum\limits_{l=0}^{i}i!\frac{S(k,i-l)}{l!}\sum\limits_{j=0}^LL!\frac{S(i,L-j)}{j!}$

　　然后我们神奇的发现后面的两个式子形式相同，以后一个为例：

　　　　$\sum\limits_{j=0}^LL!\frac{S(i,L-j)}{j!}$，换成枚举$L-j$

　　　　$=\sum\limits_{j=0}^LL!\frac{S(i,j)}{(L-j)!}$

　　　　$=\sum\limits_{j=0}^LS(i,j)C_L^jj!$，这不就是结论2吗，直接等于$L^i$

　　因此原式直接变为$\sum\limits_{i=0}^nL^ii^k$，但$n$依然很大。

　　注意当$L=1$时式子变为$\sum\limits_{i=0}^ni^k$，需要用拉格朗日插值法求解，详见CF622F.

　　若$L\neq1$，我的做法是把$i^k$变回斯特林数，即$\sum\limits_{i=0}^nL^i\sum\limits_{j=0}^kS(k,j)C_i^jj!$

　　交换求和号，即$\sum\limits_{j=0}^kS(k,j)j!\sum\limits_{i=j}^nC_i^jL^i$，考虑如何快速求出$\sum\limits_{i=j}^nC_i^jL^i$。

　　设$g(x)=\sum\limits_{i=x}^nC_i^xL^i$

　　推导一波：

　　$g(x)=\sum\limits_{i=x}^nL^iC_i^x$
　　$=\sum\limits_{i=x}^nL^iC_{i-1}^{x-1}+\sum\limits_{i=x}^nL^iC_{i-1}^x$
　　$=L\sum\limits_{i=x-1}^{n-1}L^{i}C_{i}^{x-1}+L\sum\limits_{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)$

　　到此为止，原式变为$\sum\limits_{j=1}^kS(k,j)j!g(j)$，问题就剩下如何求第二类斯特林数$S(k,j)$了。

　　应用第三个结论：$S(n,m)=\frac{1}{m!}\sum\limits_{i=1}^m(-1)^iC_m^i(m-i)^n$，将其化为卷积形式：

　　$S(n,m)=\sum\limits_{i=1}^m\frac{(-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;
}