Fibonacci Sum 推公式

$F(n)$为Fibonacci数列第n项，求$\sum (F(iC{)}^k,0\le i\le N$。
给多组N,C,K,$1\le N,C\le 10^{18}$,$k\le 10^5$

解法：

把斐波那契数列拆成这样
$$F_n=d(a^n-b^n)$$
$$d=\frac{1}{\sqrt{5}},a=\frac{1+\sqrt{5}}{2},b=\frac{1-\sqrt{5}}{2}$$
这三个数在题目里不会改变，先使用cipolla求二次剩余以及费马小定理求逆元都求出来，设为常量，我的代码被不设为常量就被卡超时了。
接下来:
$$\sum _{i=0}^N(F(iC{)}^k)=d^k\sum _{i=0}^N\sum _{j=0}^k{C}_k^j((a^{ic}{)}^j(-b^{ic}{)}^{k-j})$$
$$=d^k\sum _{j=0}^k(-1{)}^{k-j}{C}_k^j\sum _{i=0}^N(a^j{b}^{k-j}{)}^{ic}$$
令$X_j=(a^j{b}^{k-j}{)}^c$,原式变为：
$$d^k\sum _{j=0}^k{C}_k^j(-1{)}^{k-j}\sum _{i=0}^N{X}_j^i$$
并且，有等比数列求和：$$\sum _{i=0}^N{X}_j^i=\frac{X_j^{N+1}-1}{X_j-1}$$那么最终，$$d^k\sum _{j=0}^k{C}_k^j(-1{)}^{k-j}\frac{X_j^{N+1}-1}{X_j-1}$$
把阶乘，组合数都预处理一下，需要N,C,K的放在里面计算，N和C特别大，对模数欧拉降幂（就是快速幂的时候my_pow(a,b%(p-1)）
时间复杂度大概是$200\ast O(klogk)$

hduoj里提交不能在程序里面跑二次剩余以及多余的求逆元操作，貌似时间卡得很死。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <cstdlib>
using namespace std;
typedef long long ll;

const int N = 1e5+5;
const ll mod = 1e9+9;
const ll SQRT = 383008016, RSQRT = 276601605;
const ll A = 691504013, RA = 691504012;
const ll B = 308495997, RB = 308495996;
ll fac[N], in[N];
ll n, c, k;


inline ll ksm(ll a, ll b){
	ll ans = 1%mod;
	while(b){
		if(b&1) ans = ans*a%mod;
		a = a*a%mod;
		b >>= 1;
	}
	return ans%mod;
}

inline ll inv(ll x){
	return ksm(x, mod-2);
}

ll C(ll a, ll b){
	if(a<0 || b > a) return 0;
	return fac[a]*in[b]%mod*in[a-b]%mod;
} 

void init(int n){
	fac[0] = 1;
	for(int i = 1; i <= n; i++) fac[i] = fac[i-1]*(ll)i%mod;
	in[n] = inv(fac[n]);
	for(int i = n-1; i >= 0; i--) in[i] = (ll)(i+1)*in[i+1]%mod;
}
int main(){
	init(1e5);
	int _;
	scanf("%d", &_);
	while(_--){
		scanf("%lld%lld%lld", &n, &c, &k);
		ll sum = 0;
		ll bc = ksm(B, c%(mod-1)), ac = ksm(A, c%(mod-1));
		ll inv_bc = inv(bc);
		ll base = ksm(bc, k), mul = ac*inv_bc%mod;
		for(int i = 0; i <= k; i++){
			ll tmp1, tmp2;
			if(base == 1){
				tmp1 = n%mod;
			}
			else tmp1 = base*(ksm(base, n%(mod-1))+mod-1)%mod*inv(base-1)%mod;
			tmp2 = C(k, i)*tmp1%mod;
			if((k-i)&1){
				sum -= tmp2;
				sum = (sum+mod)%mod;
			}
			else{
				sum = (sum+tmp2)%mod;
			}
			base = base*mul%mod;
		}
		sum = sum*ksm(RSQRT, k)%mod;
		printf("%lld\n", sum);
	}
	return 0;
} 
/*
2
5 1 1
2 1 2
*/