【Fibonacci 序列+第一类Stirling数+二项式定理】[CodeForces - 717A]Festival Organization

题目大意

问有多少种方案选择k个不同的长度相同01串。
这些01串中要求不能出现连续的两个0。长度在[l,r]区间内。

分析

很容易发现，长度为i合法01串个数为Fi+2(Fi表示斐波那契数列的第i项)，方案数就为(Fi+2k)，令Sn=∑n+2i=3(Fik)，则答案就是Sr−Sl−1
我们来看一下怎么计算Sn

Sn=∑i=3n+2(Fik)=∑i+3n+2Fi!k!(Fi−k)!=1k!∑i=3n+2Fi(Fi−k)!=1k!∑i=3n+2∏j=0k−1Fi−j=1k!∑i=3n+2Fk−i

xk−是a的下降幂，等于∏k−1i=0x−i，也可以用有符号第一类斯特林数来化简。
xk−=∑ki=0S(k,i)xk

我们就可以继续化简

Sn=1k!∑i=3n+2∑j=0kS(k,j)Fij=1k!∑j=0kS(k,j)∑i=3n+2Fji

我们只需要快速求出∑n+2i=3Fji，矩阵好像并不能求。
考虑斐波那契数列的通项公式。

Fi=15√[(1+5√2)i−(1−5√2)i]

令ai=(1+5√2)i bi=(1−5√2)i

Fi=15√[ai−bi]Fij=15√[ai−bi]j

当我们发现将Fij用二项式定理展开之后，我们可以对同类项分别求和，可以发现这是一个等比数列求和，可以O(logn)（快速幂）计算出。
枚举{F_i}的次数，再枚举同类项，再求和，时间复杂度为O(k2logn)。
通项公式带根号无法取模，在计算的时候把数字表示成a+b5√即可。

代码

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
#define MOD 1000000007
#define MAXK 200
int k,C[MAXK+10][MAXK+10],S[MAXK+10][MAXK+10],ans;
LL l,r;
void prepare(){
    int i,j;
    for(i=0;i<=k;i++){
        S[i][i]=1;
        C[i][0]=1;
        for(j=1;j<=i;j++){
            C[i][j]=(C[i-1][j-1]+C[i-1][j])%MOD;
            S[i][j]=((S[i-1][j-1]-1ll*S[i-1][j]*(i-1))%MOD+MOD)%MOD;
        }
    }
}
template<class T>
void Read(T &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
inline void read(){
    Read(k),Read(l),Read(r);
    //f[i][0]=f[i-1][1]
    //f[i][1]=f[i-1][0]+f[i][1]
    //f[i][0]+f[i][1]=F[i+2]
    //F[i]=(((1+sqrt(5))/2)^n-((1-sqrt(5))/2)^n)/sqrt(5)
}
int quick_pow(int a,int b){
    int ret(1);
    while(b){
        if(b&1)
            ret=1ll*ret*a%MOD;
        a=1ll*a*a%MOD;
        b>>=1;
    }
    return ret;
}
struct data{
    int a,b; //a+bsqrt(5)
    inline data(){
    }
    inline data(int a,int b=0):a(a),b(b){
    }
    inline data operator + (const data &x)const{
        return data((a+x.a)%MOD,(b+x.b)%MOD);
    }
    inline data & operator +=(const data &x){
        return *this=*this+x;
    }
    inline data operator - (const data &x)const{
        return data((a-x.a+MOD)%MOD,(b-x.b+MOD)%MOD);
    }
    inline data & operator -=(const data &x){
        return *this=*this-x;
    }
    inline data operator * (const data &x)const{
        return data((1ll*a*x.a+5ll*b*x.b)%MOD,(1ll*a*x.b+1ll*b*x.a)%MOD);
    }
    inline data & operator *= (const data &x){
        return *this=*this*x;
    }
    inline data operator / (const data &x)const{
        static data ret;
        int inv;
        ret=*this*data(x.a,MOD-x.b);
        inv=quick_pow(((1ll*x.a*x.a-5ll*x.b*x.b)%MOD+MOD)%MOD,MOD-2);
        ret.a=1ll*ret.a*inv%MOD;
        ret.b=1ll*ret.b*inv%MOD;
        return ret;
    }
    inline data & operator /= (const data &x){
        return *this=*this/x;
    }
};
data quick_pow(data a,LL b){
    data ret(1);
    while(b){
        if(b&1)
            ret*=a;
        a*=a;
        b>>=1;
    }
    return ret;
}
int cal(LL n,int m){
    int i;
    data ret(0);
    data a((MOD+1)>>1,(MOD+1)>>1),b((MOD+1)>>1,1ll*(MOD-1)*((MOD+1)>>1)%MOD),gb,st;
    for(i=0;i<=m;i++){
        gb=quick_pow(a,i)*quick_pow(b,m-i);
        st=quick_pow(gb,3);
        data t=(quick_pow(gb,n+3)-st)/(gb-1);
        if(gb.a==1&&gb.b==0){
            t=n%MOD;
        }
        if((m-i)&1)
            ret-=t*C[m][i];
        else
            ret+=t*C[m][i];
    }
    ret/=quick_pow(data(0,1),m);
    return ret.a;
}
void solve(){
    int i,fac=1;
    for(i=1;i<=k;i++)
        ans=((ans+1ll*S[k][i]*(cal(r,i)-cal(l-1,i)))%MOD+MOD)%MOD;
    ans=(ans+(r-l+1)%MOD*S[k][0])%MOD;
    for(i=1;i<=k;i++)
        fac=1ll*i*fac%MOD;
    ans=1ll*ans*quick_pow(fac,MOD-2)%MOD;
}
int main()
{
    read();
    prepare();
    solve();
    printf("%d\n",ans);
}