codeforces 717A. Festival Organization

题目链接：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;
}