Fibonacci Sum
题目链接
http://acm.hdu.edu.cn/showproblem.php?pid=6755
思路
首先关于题目需要的知识点单独记录在了下面的博客
https://blog.youkuaiyun.com/xmyrzb/article/details/107526823
求出斐波那契数列的通项公式为 F n = d ( a n − b n ) F_n = d(a^n-b^n) Fn=d(an−bn),其中 d = 1 5 \quad d=\frac 1 {\sqrt 5} d=51, a = 1 + 5 2 \quad a=\frac{1+\sqrt 5}{2} a=21+5, b = 1 − 5 2 \quad b=\frac{1-\sqrt 5}{2} b=21−5,在这里可以提前算出这些值(用逆元和二次剩余计算)
关于sq5用二次剩余算出来是 383008016 383008016 383008016, i n v inv inv表示计算逆元
a = 691504013 , b = 308495997 a=691504013,b=308495997 a=691504013,b=308495997,
a=(1+sq5)*inv(2)%mod;
b=(1-sq5+mod)%mod*inv(2)%mod;
可以算出通项公式为
d k ∑ j = 0 k C k j x n + 1 − 1 x − 1 , x = ( a j b k − j ) C d^k \sum_{j=0}^kC_k^j \frac {x^{n+1}-1}{x-1}, \quad x=(a^jb^{k-j})^C dk∑j=0kCkjx−1xn+1−1,x=(ajbk−j)C
然后需要求幂的都使用欧拉降幂,指数对1e9+8取模
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1e5+5;
const ll sqrt5=383008016;
const ll A=691504013,B=308495997;
const ll P=1e9+9;
ll fac[maxn],finv[maxn];
ll fp(ll b,ll p){ //快速幂
ll ans=1;
while(p){
if(p&1){
ans=ans*b%P;
}
p>>=1;
b=b*b%P;
}
return ans;
}
void init(int n){
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%P; //计算阶乘
finv[n]=fp(fac[n],P-2); //计算阶乘的逆元
for(int i=n-1;i>=0;i--) finv[i]=finv[i+1]*(i+1)%P;
// finv[1] = 1; //计算阶乘逆元的第二种方法
// for( int i = 2; i <=n; i++ ) finv[i] = ( P - P / i ) * finv[ P % i ] % P;
// finv[0]=1;
// for(int i=1;i<=n;i++) finv[i]=finv[i]*finv[i-1]%P;
}
ll C(ll n,ll m){ //计算系数C(k,i)
if(m<0 || m>n) return 0;
return fac[n]*finv[n-m]%P*finv[m]%P; // !n/!(n-m) / !m
}
ll solve(ll n,ll c,ll k){
ll A_C=fp(A,c%(P-1)); //用欧拉降幂公式和快速幂求解A^C, B^C
ll B_C=fp(B,c%(P-1));
ll inv_B_C=fp(B_C,P-2); //求B^C的逆元(使用费马小定律,快速幂
ll base=fp(B_C,k); //求得初始时的x,x为A^0 * B^k
ll multi=A_C*inv_B_C%P; //求得每次循环时x变化的量
ll ans=0; //最终值
for(ll i=0;i<=k;i++){
ll temp;
//首先计算(x^n-1)/(x-1)
if(base==1){ //特判一下x的值,若为1则值为n,为1的时候分母会为0
temp=n%P;
}
else{ //求(x^n-1)/(x-1),欧拉降幂求x^n,费马小定律求(x-1)的逆元
temp=base*(fp(base,n%(P-1))-1+P)%P*fp(base-1,P-2)%P;
}
//将系数C(k,j)乘上
ll temp2=C(k,i)*temp%P;
//判断正负性
if((k-i)&1){
ans-=temp2;
ans=(ans+P)%P; //保证值在P之内
}
else{
ans+=temp2;
ans=(ans+P)%P;
}
//将值乘上每次循环的变量,得到下一次的x
base=base*multi%P;
}
ll inv_sqrt5=fp(sqrt5,P-2); //求根号5的逆元
ans=ans*fp(inv_sqrt5,k)%P; //将根号5分之一的k次方乘上去
return ans;
}
int main () {
init(1e5); //提前对阶乘和阶乘逆元打表
int T;
scanf("%d",&T);
while(T--){
ll n,c,k;
scanf("%lld%lld%lld",&n,&c,&k);
ll ans=solve(n,c,k);
printf("%lld\n",ans);
}
}
代码取自:
作者: UCPRER
出处:https://www.cnblogs.com/ucprer/p/13361712.html