51nod1236 序列求和V3

题目链接

题目

给定$n,k$，求$S(n,k)={\sum }_{i=1}^n{f}_i^k\mathrm{m}\mathrm{o}\mathrm{d}1000000009$，其中$f_i$是斐波那契数列的第$i$项。
$(1\le n\le 10^{18},1\le k\le 10^5)$

题解：

${\sum }_{i=1}^n{f}_i^k={\sum }_{i=1}^n{\left(A(B^i-C^i)\right)}^k$，其中$A=\frac{\sqrt{5}}{5},B=\frac{1+\sqrt{5}}{2},C=\frac{1-\sqrt{5}}{2}$。

$S(n,k)=A^k{\sum }_{i=1}^n(B^i-C^i{)}^k$

$=A^k{\sum }_{i=1}^n{\sum }_{j=0}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right)(-1{)}^{k-j}{B}^{ik}{C}^{i(k-j)}$

$=A^k{\sum }_{j=0}^k(-1{)}^{k-j}(\genfrac{}{}{0ex}{}{k}{j}){\sum }_{i=1}^n(B^k{C}^{k-j}{)}^i$

后半部分就是一个等比数列求和，所以$O(k)$计算即可。注意公比为1的等比数列不可用求和公式计算，$\sqrt{5}$在模1000000009意义下的数用二次剩余计算即可。

复杂度：$O(klogn)$

代码：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<vector>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<string>
#include<bitset>
#include<sstream>
#include<ctime>
//#include<chrono>
//#include<random>
//#include<unordered_map>
using namespace std;

#define ll long long
#define ls o<<1
#define rs o<<1|1
#define pii pair<int,int>
#define fi first
#define se second
#define pb push_back
#define mp make_pair
#define sz(x) (int)(x).size()
#define all(x) (x).begin(),(x).end()
const double pi=acos(-1.0);
const double eps=1e-6;
const int mod=1e9+9;
const int INF=0x3f3f3f3f;
const int maxn=1e5+5;
ll read(){
	ll x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
ll qpow(ll a,ll p=mod-2){
	ll res=1;
	while(p){
		if(p&1)res=res*a%mod;
		a=a*a%mod;
		p>>=1;
	}
	return res;
}
ll n;
ll sqrt5,A,B,C;
int k;
ll fac[maxn],inv[maxn],Bn[maxn],Cn[maxn];
ll fix(ll x){
	return (x%mod+mod)%mod;
}
void init(int n){
	sqrt5=383008016;
	A=qpow(sqrt5)%mod;
	B=(1+sqrt5)*qpow(2)%mod;
	C=fix((1-sqrt5)*qpow(2));
	fac[0]=1;    
    for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;   
    inv[n]=qpow(fac[n]);    
    for(int i=n-1;i>=0;i--){    
        inv[i]=inv[i+1]*(i+1)%mod;    
    }    
    Bn[0]=Cn[0]=1;
    for(int i=1;i<=n;i++){
    	Bn[i]=Bn[i-1]*B%mod;
    	Cn[i]=Cn[i-1]*C%mod;
    }
}   
ll com(int n,int m){    
    if(n<m)return 0;  
    return fac[n]*inv[m]%mod*inv[n-m]%mod;    
}    
ll S(ll x,ll n){
	if(x==1)return n%mod;
	else return x*(qpow(x,n%(mod-1))-1)%mod*qpow(x-1)%mod;
}
int main(void){
	// freopen("in.txt","r",stdin);
	init(100000);
	int T;
	scanf("%d",&T);
	while(T--){
		scanf("%lld%d",&n,&k);
		ll Ak=qpow(A,k);
		ll ans=0,f=1;
		for(int i=k;i>=0;i--){
			ll tmp=Bn[i]*Cn[k-i]%mod;
			ll cur=S(tmp,n);
			cur=cur*com(k,i)%mod;
			ans=(ans+f*cur)%mod;
			f*=-1;
		}
		ans=(ans+mod)%mod;
		ans=ans*Ak%mod;
		printf("%lld\n",ans);
	}
	return 0;
}