HDOJ 6755 - Fibonacci Sum

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

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

\(n,c\in\left[1,10^{18}\right],k\in\left[1,10^5\right],T\in[1,200]\)。\(\mathrm{TL}=3\mathrm s\)。

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

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

带到原式里面得

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

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

\[ans=\left(\frac1{\sqrt5}\right)^k\sum_{i=0}^n\left(a^{ic}-b^{ic}\right)^k \]

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

二项式定理展开得

\[\begin{aligned}ans&=\left(\frac1{\sqrt5}\right)^k\sum_{i=0}^n\sum_{j=0}^k\binom kja^{icj}(-b^{ic})^{k-j}\\&=\left(\frac1{\sqrt5}\right)^k\sum_{i=0}^n\sum_{j=0}^k(-1)^{k-j}\binom kja^{icj}b^{ic(k-j)}\end{aligned} \]

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

\[ans=\left(\frac1{\sqrt5}\right)^k\sum_{i=0}^k(-1)^{k-i}\binom ki\sum_{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(Tk\log 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<<"\n";
}
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;
}

UPD 2020.8.27：偶然看到队爷 lqs 的博客里这题的题解，发现可以线性求一堆数的逆元？

然后到洛谷上翻了模板，发现我今年 1.13 就 A 掉了，但是现在却一点关于线性求逆元的印象都没有？可见当时我是多么不求甚解。

进入正题。这个东西其实很简单……要求 \(n\) 个数 \(a_i\) 的逆元的话，考虑将所有数的乘积的逆元 \(\Pi=\prod a_i^{-1}\) 给算出来，那么 \(a_i^{-1}=\Pi\cdot\prod\limits_{j=1}^{i-1}a_j\cdot\prod\limits_{j=i+1}^na_j\)，这样一来把它变成乘法的事情就好办多了，这两个 \(\prod\) 只需要前后缀积即可线性求出。显然只求了一次逆元，其他操作都是线性，那么总时间复杂度 \(\mathrm O(\log p+n)\)。

回到本题。只需要将所有希望求逆元的数给离线下来一起求，即可线性。复杂度看实现，最好的是 \(\mathrm O(\log p+Tk)\)，为了方便我写的是 \(\mathrm O(T(\log p+k))\)。

代码（\(920\mathrm{ms}\to 234\mathrm{ms}\)）：

#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;}
int prod_inv,Prod[K+2],proD[K+3];
int _inv(int x){return prod_inv*Prod[x]%mod*proD[x+2]%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,
		frac_cpy=frac;
	for(int i=0;i<=k;i++){
		if(frac-1)Prod[i+1]=frac-1;
		else Prod[i+1]=1;
		(((frac*=frac_mul)%=mod)+=mod)%=mod; 
	}
	frac=frac_cpy;
	Prod[0]=proD[k+2]=1;
	for(int i=k+1;i;i--)proD[i]=proD[i+1]*Prod[i]%mod;
	for(int i=1;i<=k+1;i++)Prod[i]=Prod[i]*Prod[i-1]%mod;
	prod_inv=inv(Prod[k+1]);
	for(int i=0;i<=k;i++){
		int s=(frac==1?n+1:(ed*frac%mod-1)*_inv(i))%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<<"\n";
}
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;
}