题目链接：CF717A
翻译：luogu
对于一个确定的长度(n)，合法的方案数为(fib_{n+2})
所以最后求的就是(sum_{i=l}^rdbinom{fib_{i+2}}{k})
记(f_n=sum_{i=0}^ndbinom{fib_i}{k})，那么答案也就是(f_{r+2}-f_{l+1})
推一下式子
[egin{aligned}
f_n=&sum_{i=0}^nfrac{fib_i!}{k!(fib_i-k)!}\
=&frac{1}{k!}sum_{i=0}^nfrac{fib_i!}{(fib_i-k)!}\
=&frac{1}{k!}sum_{i=0}^nfib_i^{underline{k}}\
=&frac{1}{k!}sum_{i=0}^nsum_{j=0}^k(-1)^{j-k}egin{bmatrix}k\jend{bmatrix}fib_i^j\
=&frac{1}{k!}sum_{j=0}^kegin{bmatrix}k\jend{bmatrix}(-1)^{k-j}sum_{i=0}^nfib_i^j
end{aligned}
]
基本上就是运用了一下第一类斯特林数和下降幂之间的关系
问题就是对每个给定的(j)，怎么求(sum_{i=0}^n fib_i^j)
我们有(fib_n=frac{(frac{1+sqrt5}{2})^n-(frac{1-sqrt5}{2})^n}{sqrt5})，令(a=frac{1+sqrt5}{2},b=frac{1-sqrt5}{2})
则(sum_{i=0}^nfib_i^j=sum_{i=0}^nfrac{(a^i-b^i)^j}{(sqrt5)^j}=frac{sum_{k=0}^j {jchoose k}(-1)^ksum_{i=0}^n(b^ka^{j-k})^i}{(sqrt5)^j})
上面这个式子用二项式定理展开后，最后的那个sigma可以用等比数列求和解决
貌似时间是(O(klogr))，结合一开始的最外面枚举(j)总时间是(O(k^2logr))？
时间没有问题，但是注意到(5)在模(1000000007)意义下并没有二次剩余
考虑扩域，即类比虚数，将所有数定义成(a+bsqrt5)的形式，重新定义四则运算即可（除法的话直接考虑分母有理化即可）
#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
#include<queue>
#include<map>
#include<set>
using namespace std;
#define lowbit(x) (x)&(-x)
#define sqr(x) (x)*(x)
#define fir first
#define sec second
#define rep(i,a,b) for (register int i=a;i<=b;i++)
#define per(i,a,b) for (register int i=a;i>=b;i--)
#define maxd 1000000007
#define inv2 500000004
#define eps 1e-6
typedef long long ll;
const int N=200;
const double pi=acos(-1.0);
int k;
ll l,r,c[220][220],s[220][220],invk=1;
ll qpow(ll x,ll y)
{
ll ans=1;
while (y)
{
if (y&1) ans=ans*x%maxd;
x=x*x%maxd;y>>=1;
}
return ans;
}
struct Complex{
ll x,y;
Complex(ll _x=0,ll _y=0) {x=_x;y=_y;}
};
Complex operator +(Complex a,Complex b)
{
return Complex((a.x+b.x)%maxd,(a.y+b.y)%maxd);
}
Complex operator -(Complex a,Complex b)
{
return Complex((a.x-b.x+maxd)%maxd,(a.y-b.y+maxd)%maxd);
}
Complex operator *(Complex a,Complex b)
{
return Complex((a.x*b.x+a.y*b.y%maxd*5)%maxd,(a.x*b.y+a.y*b.x)%maxd);
}
Complex operator /(Complex a,Complex b)
{
a=a*Complex(b.x,maxd-b.y);b=b*Complex(b.x,maxd-b.y);
ll inv=qpow(b.x,maxd-2);
a.x=a.x*inv%maxd;a.y=a.y*inv%maxd;
return a;
}
Complex qpow(Complex a,ll y)
{
Complex ans=Complex(1,0);
while (y)
{
if (y&1) ans=ans*a;
a=a*a;y>>=1;
}
return ans;
}
void init()
{
s[0][0]=1;c[0][0]=1;
rep(i,1,N)
{
c[i][0]=1;
rep(j,1,i)
{
c[i][j]=(c[i-1][j-1]+c[i-1][j])%maxd;
s[i][j]=(s[i-1][j-1]+s[i-1][j]*(i-1))%maxd;
}
}
rep(i,1,k) invk=invk*i%maxd;
invk=qpow(invk,maxd-2);
}
Complex calc(Complex a,ll n)
{
Complex tmp=Complex(1,0);
if((a.x==1) && (!a.y)) return Complex((n+1)%maxd,0);
else return (qpow(a,n+1)-tmp)/(a-tmp);
}
ll solve(ll n)
{
Complex a=Complex(inv2,inv2),b=Complex(inv2,maxd-inv2);
ll ans=0;
rep(i,0,k)
{
ll now=s[k][i];
if ((k-i)&1) now=(maxd-now)%maxd;
Complex sum=Complex(0,0);
rep(j,0,i)
{
Complex tmp1=Complex(0,0);
if (j&1) tmp1.x=maxd-c[i][j];else tmp1.x=c[i][j];
tmp1=tmp1*calc(qpow(b,j)*qpow(a,i-j),n);
sum=sum+tmp1;
}
Complex inv=Complex();
if (i&1) inv.y=qpow(5,i/2);else inv.x=qpow(5,i/2);
sum=sum/inv;
ans=(ans+now*sum.x)%maxd;
}
ans=ans*invk%maxd;
return ans;
}
int main()
{
scanf("%d%lld%lld",&k,&l,&r);
init();
printf("%lld",(solve(r+2)-solve(l+1)+maxd)%maxd);
return 0;
}