CF717A Festival Organization（第一类斯特林数，斐波那契数列）

Posted 2022-01-15 1000suns

题目大意：求 $\sum\limits_n=l^r\dbinomf_nk\bmod 10^9+7$。其中 $f_n$ 是长度为 $n$ 的 $01$ 序列中，没有连续两个或超过两个 $0$ 的个数。

$1\le k\le 200,1\le l\le r\le 10^18$。

---

先考虑如何求 $f_n$。

令 $g[i][j]$ 表示长度为 $i$，结尾是 $j$ 的序列个数。

$$g[i][0]=g[i-1][1]$$

$$g[i][1]=g[i-1][0]+g[i-1][1]$$

将第一个式子代入第二个式子有 $g[i][1]=g[i-2][1]+g[i-1][1]$。

手玩发现 $g[1][1]=1,g[2][1]=2$。那么有 $g[i][1]=fib_i+1$。（$fib$ 是斐波那契数列）

要求 $g[i][0]+g[i][1]=g[i-1][1]+g[i][1]=g[i+1][1]=fib_i+2$。

那么式子就是 $\sum\limits_n=l+2^r+2\dbinomfib_nk$。（为方便 $l+=2,r+=2$，下文假设是 $l$ 到 $r$）

$$\sum\limits_n=l^r\dbinomfib_nk$$

$$\dfrac1k!\sum\limits_n=l^rfib_n^\underlinek$$

$$\dfrac1k!\sum\limits_n=l^r\sum\limits_i=0^k(-1)^k-i\beginbmatrixk\\i\endbmatrixfib_n^i$$

$$\dfrac1k!\sum\limits_i=0^k(-1)^k-i\beginbmatrixk\\i\endbmatrix\sum\limits_n=l^rfib_n^i$$

前面可以随便枚举，然而斐波那契 $k$ 次方前缀和这东西怎么搞？

这时就要用上大名鼎鼎的通项公式：（没记住的也可以现场用特征方程推）

$$fib_n=\dfrac\sqrt55(\dfrac1+\sqrt52)^n-\dfrac\sqrt55(\dfrac1-\sqrt52)^n$$

令 $a=\dfrac\sqrt55,b=-\dfrac\sqrt55,x=\dfrac1+\sqrt52,y=\dfrac1-\sqrt52$，那么 $fib_n=ax^n+by^n$。

$$\dfrac1k!\sum\limits_i=0^k(-1)^k-i\beginbmatrixk\\i\endbmatrix\sum\limits_n=l^r(ax^n+by^n)^i$$

$$\dfrac1k!\sum\limits_i=0^k(-1)^k-i\beginbmatrixk\\i\endbmatrix\sum\limits_n=l^r\sum\limits_j=0^i\dbinomij(ax^n)^j(by^n)^i-j$$

$$\dfrac1k!\sum\limits_i=0^k(-1)^k-i\beginbmatrixk\\i\endbmatrix\sum\limits_n=l^r\sum\limits_j=0^i\dbinomij(a^jb^i-j)(x^jy^i-j)^n$$

$$\dfrac1k!\sum\limits_i=0^k(-1)^k-i\beginbmatrixk\\i\endbmatrix\sum\limits_j=0^i\dbinomij(a^jb^i-j)\sum\limits_n=l^r(x^jy^i-j)^n$$

此时 $i$ 和 $j$ 可以 $O(k^2)$ 枚举，而 $n$ 那里是个等差数列，总复杂度 $O(k^2\log r)$。（可能要预处理出所有斯特林数和组合数才能做到这个复杂度）

但是有个严重的问题：$\sqrt5$ 在模 $10^9+7$ 意义下不存在。

那么就要用到一个骚操作：扩系。

令 $(a,b)=a+b\sqrt5$，那么上文中的四个常数 $a=(0,5^-1),b=(0,-5^-1),x=(2^-1,2^-1),y=(2^-1,-2^-1)$。

那么进行加减乘除就简单了：

$$(a,b)\pm(c,d)=(a\pm c,b\pm d)$$

$$(a,b)\times(c,d)=(ac+5bd,ad+bc)$$

$$\dfrac1(a,b)=(\dfracaa^2-5b^2,\dfrac-ba^2-5b^2)$$

因为式子没推错（只能这么想了啊），最后求出的一定是整数。

那么就做完了。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=222,mod=1000000007,inv2=500000004,inv5=400000003;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline ll read()
    char ch=getchar();ll x=0,f=0;
    while(ch<‘0‘ || ch>‘9‘) f|=ch==‘-‘,ch=getchar();
    while(ch>=‘0‘ && ch<=‘9‘) x=x*10+ch-‘0‘,ch=getchar();
    return f?-x:x;

int k,fac[maxn],invfac[maxn],ans,S[maxn][maxn],C[maxn][maxn];
ll l,r;
inline int add(int x,int y)return x+y<mod?x+y:x+y-mod;
inline int sub(int x,int y)return x<y?x-y+mod:x-y;
inline int mul(int x,int y)return 1ll*x*y%mod;
inline int qpow(int a,ll b)
    int ans=1;
    for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a);
    return ans;

struct comp
    int x,y;
    comp(const int xx=0,const int yy=0):x(xx),y(yy)
    inline comp operator+(const comp &c)constreturn comp(add(x,c.x),add(y,c.y));
    inline comp operator-(const comp &c)constreturn comp(sub(x,c.x),sub(y,c.y));
    inline comp operator*(const comp &c)constreturn comp(add(mul(x,c.x),mul(5,mul(y,c.y))),add(mul(x,c.y),mul(y,c.x)));
    inline comp inv()const
        comp ans(x,y?mod-y:0);
        int dn=qpow(sub(mul(x,x),mul(5,mul(y,y))),mod-2);
        return ans*dn;
    
    inline comp operator/(const comp &c)constreturn *this*c.inv();
    inline bool operator==(const comp &c)constreturn x==c.x && y==c.y;
a(0,inv5),b(0,mod-inv5),x(inv2,inv2),y(inv2,mod-inv2);
inline comp cqpow(comp a,ll b)
    comp ans(1,0);
    for(;b;b>>=1,a=a*a) if(b&1) ans=ans*a;
    return ans;

comp calc(comp x,ll l,ll r)
    if(x==1) return (r-l+1)%mod;
    return (cqpow(x,r+1)-cqpow(x,l))/(x-1);

int main()
    k=read();l=read();r=read();
    FOR(i,0,k) C[i][0]=C[i][i]=1;
    FOR(i,1,k) FOR(j,1,i-1) C[i][j]=add(C[i-1][j],C[i-1][j-1]);
    S[1][1]=1;
    FOR(i,2,k) FOR(j,1,i) S[i][j]=add(mul(i-1,S[i-1][j]),S[i-1][j-1]);
    fac[0]=1;
    FOR(i,1,k) fac[i]=mul(fac[i-1],i);
    invfac[k]=qpow(fac[k],mod-2);
    FOR(i,0,k)
        int s=0;
        FOR(j,0,i)
            comp tmp1=cqpow(a,j)*cqpow(b,i-j),tmp2=cqpow(x,j)*cqpow(y,i-j);
            s=add(s,mul(C[i][j],(tmp1*calc(tmp2,l+2,r+2)).x));
        
        s=mul(s,S[k][i]);
        if((k-i)&1) ans=sub(ans,s);
        else ans=add(ans,s);
    
    printf("%d\n",mul(ans,invfac[k]));

View Code