牛客多校第九场 && ZOJ3774 The power of Fibonacci（二次剩余定理+斐波那契数列通项/循环节）题解
题意1.1：
求\(\sum_{i=1}^n Fib^m\mod 1e9+9\)，\(n\in[1, 1e9], m\in[1, 1e4]\)
思路1.1
我们首先需要知道斐波那契数列的通项是：\(Fib_i = \frac{\sqrt5}{5}[(\frac{1+\sqrt5}{2})^i-(\frac{1-\sqrt5}{2})^i]\)，因为取模是个质数，我们可以用二次剩余定理得到\(\sqrt5 \mod 1e9+9 = 383008016\)，然后就可以得到\(\frac{\sqrt5}{5}, \frac{1+\sqrt5}{2},\frac{1-\sqrt5}{2}\)的取模的整数值，我们记为\(s =\frac{\sqrt5}{5}, r_1=\ \frac{1+\sqrt5}{2},r_2= \frac{1-\sqrt5}{2}\)。那么
\[\sum_{i=1}^nFib_i^m = \sum_{i=1}^ns^m(r_1^i - r_2^i)^m = s^m\sum_{i=1}^n\sum_{r=0}^m[(-1)^rC_m^rr_1^{i(m-r)}r_2^{ir}]
\]
因为\(n\)太大了，所以我们还要继续化简：
\[ s^m\sum_{i=1}^n\sum_{r=0}^m[(-1)^rC_m^rr_1^{i(m-r)}r_2^{ir}] = s^m\sum_{r=0}^m[(-1)^rC_m^r\sum_{i=1}^n(r_1^{i(m-r)}r_2^{ir})]
\]
因为\(\sum_{i=1}^n(r_1^{i(m-r)}r_2^{ir})\)是个等比数列，故我们假设\(q_r = r_1^{m-r}r_2^{r}\)，则由等比数列性质可得：
\[s^m\sum_{r=0}^m[(-1)^rC_m^r\sum_{i=1}^n(r_1^{i(m-r)}r_2^{ir})]= s^m\sum_{r=0}^m[(-1)^rC_m(q_r\frac{1-q_r^n}{1-q_r})]
\]
推到这里就结束了，直接求解。
题意1.2
牛客传送门
求\(\sum_{i=1}^n Fib^m\mod 1e9\)，\(n\in[1, 1e9], m\in[1, 1e4]\)
思路1.2
因为取模是个合数，那就不能二次剩余定理了。斐波那契数列的取模是有循环节的，斐波那契数列幂次的取模循环节和原数列取模循环节相同，那么我们直接暴力找到循环节，然后求解，最后用中国剩余定理合在一起即可。
代码：
//zoj3774
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<ctime>
#include<cmath>
#include<cstdio>
#include<string>
#include<vector>
#include<cstring>
#include<sstream>
#include<iostream>
#include<algorithm>
#include<unordered_map>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int maxn = 100000 + 5;
const int INF = 0x3f3f3f3f;
const ull seed = 131;
const ll MOD = 1000000009;
using namespace std;
ll fac[maxn], inv[maxn];
ll ppow(ll a, ll b){
ll ret = 1;
while(b){
if(b & 1) ret = ret * a % MOD;
b >>= 1;
a = a * a % MOD;
}
return ret;
}
ll C(int n, int m){
return fac[n] * inv[m] % MOD * inv[n - m] % MOD;
}
void init(int n){
fac[0] = inv[0] = 1;
for(int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i % MOD;
inv[n] = ppow(fac[n], MOD - 2);
for(int i = n - 1; i >= 1; i--) inv[i] = (i + 1LL) * inv[i + 1] % MOD;
}
int main(){
init(100000);
ll n, m;
ll s = 276601605, r1 = 691504013, r2 = 308495997;
// printf("%lld\n", s * (r1 - r2) % MOD);
int T;
scanf("%d", &T);
while(T--){
scanf("%lld%lld", &n, &m);
ll ans = 0;
for(int r = 0; r <= m; r++){
ll q = 1LL * ppow(r1, m - r) * ppow(r2, r) % MOD;
ll sum = 1LL * q * ((ppow(q, n) - 1LL)) % MOD * ppow(q - 1, MOD - 2) % MOD;
if(q == 1) sum = n % MOD; //!!!!!
sum = 1LL * sum * C(m, r) % MOD;
if(r & 1) sum = -sum;
ans = (ans + sum) % MOD;
}
ans = (ans % MOD + MOD) % MOD;
ans = 1LL * ans * ppow(s, m) % MOD;
printf("%lld\n", ans);
}
return 0;
}
//牛客
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<ctime>
#include<cmath>
#include<cstdio>
#include<string>
#include<vector>
#include<cstring>
#include<sstream>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int maxn = 100000 + 5;
const int INF = 0x3f3f3f3f;
const ull seed = 131;
const ll MOD = 1e9;
using namespace std;
ll mod[2] = {512, 1953125}, a[2];
ll sum[7812500 + 5], fab[7812500 + 5];
int lp[2] = {768, 7812500}; // 循环节
ll ppow(ll a, ll b, ll p){
ll ret = 1;
while(b){
if(b & 1) ret = ret * a % p;
a = a * a % p;
b >>= 1;
}
return ret;
}
void exgcd(ll a, ll b, ll &x, ll &y){
if(b == 0){
x = 1, y = 0;
return;
}
exgcd(b, a % b, x, y);
ll tp = x;
x = y, y = tp - a / b * y;
}
ll CRT(){
ll ans = 0, lcm = 1, x, y;
for(int i = 0; i < 2; i++) lcm *= mod[i];
for(int i = 0; i < 2; i++){
ll tp = lcm / mod[i];
exgcd(tp, mod[i], x, y);
x = (x % mod[i] + mod[i]) % mod[i];
ans = (ans + tp * x * a[i]) % lcm;
}
return (ans % lcm + lcm) % lcm;
}
int main(){
ll n, m;
scanf("%lld%lld", &n, &m);
fab[0] = 0, fab[1] = 1, sum[0] = 0, sum[1] = 1;
for(int i = 2; i <= lp[1]; i++){
fab[i] = (fab[i - 1] + fab[i - 2]) % MOD;
sum[i] = (sum[i - 1] + ppow(fab[i], m, MOD)) % MOD;
}
a[0] = (sum[lp[0]] * (n / lp[0]) + sum[n % lp[0]]) % mod[0];
a[1] = (sum[lp[1]] * (n / lp[1]) + sum[n % lp[1]]) % mod[1];
printf("%lld\n", CRT());
return 0;
}
