HDOJ 6755

多校 Round 1的签到题之一。现在来补Round 1啥心态？鸽子的心态。

给定(n,c,k)，求(sumlimits_{i=0}^nf_{ic}^k)，其中(f_i)是斐波那契数列的第(i)项。答案对(p=10^9+9)取模。本题多测，共(T)组数据。

(n,cinleft[1,10^{18} ight],kinleft[1,10^5 ight],Tin[1,200])。(mathrm{TL}=3mathrm s)。

首先用上广为人知的斐波那契数列通项公式（其中(sqrt5)恰好是在(mod p)意义下有值的，暴力枚举一下几秒出答案：(383008016)）：

[f_i=frac1{sqrt 5}left(left(frac{1+sqrt5}2 ight)^i-left(frac{1-sqrt5}2 ight)^i ight) ]

带到原式里面得

[ans=sum_{i=0}^nleft(frac1{sqrt 5}left(left(frac{1+sqrt5}2 ight)^{ic}-left(frac{1-sqrt5}2 ight)^{ic} ight) ight)^k ]

常系数提到前面去，令(a=dfrac{1+sqrt5}2,b=dfrac{1-sqrt5}2)，得

[ans=left(frac1{sqrt5} ight)^ksum_{i=0}^nleft(a^{ic}-b^{ic} ight)^k ]

然后就可以愉快地推柿子了。

二项式定理展开得

[egin{aligned}ans&=left(frac1{sqrt5} ight)^ksum_{i=0}^nsum_{j=0}^kinom kja^{icj}(-b^{ic})^{k-j}\&=left(frac1{sqrt5} ight)^ksum_{i=0}^nsum_{j=0}^k(-1)^{k-j}inom kja^{icj}b^{ic(k-j)}end{aligned} ]

显然第一个(sum)是(mathrm O(n))的，没有前途，于是交换(sum)使得第一个(sum)是(mathrm O(k))的，得

[ans=left(frac1{sqrt5} ight)^ksum_{i=0}^k(-1)^{k-i}inom kisum_{j=0}^na^{icj}b^{(k-i)cj} ]

此时不难发现第二个(sum)是等比数列求和，首项(1)，末项(a^{icn}b^{(k-i)cn})，公比(a^{ic}b^{(k-i)c})，套公式即可，注意讨论公比为(1)的情况。

此时时间复杂度(mathrm O(Tk(log n+log c+log p)))，其中后面那一坨(log)是快速幂。要做好多次快速幂，常数非常大，不知道能不能过，反正当时现场sjc是使用了极限卡常才不T的，结果WA了。

考虑优化，尽量减少快速幂。注意到当(i)递增时，末项和公比分别组成的数列是等比数列，于是我们可以预处理它们的公比，每次递推乘以公比就可以得到下一个(i)时的末项和公比了。这样只剩下一个等比数列求和公式中分数线下面的那个式子要求逆元，而且无法递推。现在时间复杂度(mathrm O(Tklog p))，常数也很小，就能过了。

代码：

#include<bits/stdc++.h>
using namespace std; 
#define int long long
const int mod=1000000009,sqrt5=383008016;
int qpow(int x,int y){
	int res=1;
	while(y){
		if(y&1)res=res*x%mod;
		x=x*x%mod;
		y>>=1;
	}
	return res;
}
int inv(int x){return qpow(x,mod-2);}
const int K=100000;
int a,b;
int n,c,k;
int fact[K+1],factinv[K+1];
int comb(int x,int y){return fact[x]*factinv[y]%mod*factinv[x-y]%mod;}
void mian(){
	cin>>n>>c>>k;
	int ans=0;
	int frac=qpow(qpow(b,k),c),ed=qpow(qpow(qpow(b,k),c),n),//初始公比和末项 
	    frac_mul=qpow(a,c)*inv(qpow(b,c))%mod,ed_mul=qpow(qpow(a,c),n)*inv(qpow(qpow(b,c),n))%mod;//公比和末项分别组成的等比数列的公比 
	for(int i=0;i<=k;i++){
		int s=(frac==1?n+1:(ed*frac%mod-1)*inv(frac-1))%mod;
		(ans+=comb(k,i)*(k-i&1?-1:1)*s%mod)%=mod;//柿子 
		(((frac*=frac_mul)%=mod)+=mod)%=mod,(ed*=ed_mul)%=mod;//递推 
	}
	(ans*=qpow(inv(sqrt5),k))%=mod;
	cout<<(ans+mod)%mod<<"
";
}
signed main(){
	a=(1+sqrt5)*inv(2)%mod,b=(1-sqrt5)*inv(2)%mod;
	fact[0]=factinv[0]=1;
	for(int i=1;i<=K;i++)factinv[i]=inv(fact[i]=fact[i-1]*i%mod);//预处理
	int testnum;
	cin>>testnum;
	while(testnum--)mian();
	return 0;
}