codeforces 717A. Festival Organization

题目链接：CF717A

翻译：luogu

对于一个确定的长度\(n\)，合法的方案数为\(fib_{n+2}\)

所以最后求的就是\(\sum_{i=l}^r\dbinom{fib_{i+2}}{k}\)

记\(f_n=\sum_{i=0}^n\dbinom{fib_i}{k}\)，那么答案也就是\(f_{r+2}-f_{l+1}\)

推一下式子

\[\begin{aligned} f_n=&\sum_{i=0}^n\frac{fib_i!}{k!(fib_i-k)!}\\ =&\frac{1}{k!}\sum_{i=0}^n\frac{fib_i!}{(fib_i-k)!}\\ =&\frac{1}{k!}\sum_{i=0}^nfib_i^{\underline{k}}\\ =&\frac{1}{k!}\sum_{i=0}^n\sum_{j=0}^k(-1)^{j-k}\begin{bmatrix}k\\j\end{bmatrix}fib_i^j\\ =&\frac{1}{k!}\sum_{j=0}^k\begin{bmatrix}k\\j\end{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}^n\frac{(a^i-b^i)^j}{(\sqrt5)^j}=\frac{\sum_{k=0}^j {j\choose k}(-1)^k\sum_{i=0}^n(b^ka^{j-k})^i}{(\sqrt5)^j}\)

上面这个式子用二项式定理展开后，最后的那个sigma可以用等比数列求和解决

貌似时间是\(O(klogr)\)，结合一开始的最外面枚举\(j\)总时间是\(O(k^2logr)\)？

时间没有问题，但是注意到\(5\)在模\(1000000007\)意义下并没有二次剩余

考虑扩域，即类比虚数，将所有数定义成\(a+b\sqrt5\)的形式，重新定义四则运算即可（除法的话直接考虑分母有理化即可）

#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;
}