题意：
给定三个数字 \(n,c,k\) ，求以下式子
\(\sum_{i=0}^nF(ic)^k\%(10^9+9)\)
其中\(F(x)\)为斐波那契数列第\(x\)项。
\(1\leq n,c\leq10^{18},1\leq k \leq10^5\)
分析：
在比赛的时候我搜索了一波斐波那契数列的性质，意外发现了类似题目，ZOJ3774，只不过那个题里\(c=1\)，这里解释一下过程。
斐波那契有一个公式：
\[F(n)=\frac 1 {\sqrt5}((\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n) \]
我们先看\(c=1\)的情况，也就是每一项都取。即：
\[sum=\sum_{i=0}^nF(i)^k\%(10^9+9)\\F(n)^k=(\frac 1 {\sqrt5}((\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n))^k\\F(n)^k=\frac 1 {(\sqrt5)^k}((\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n)^k \]
我们假设\(a=(\frac{1+\sqrt5}{2}),b=(\frac{1-\sqrt5}{2})\)
那么有：
\[F(n)^k=\frac 1 {(\sqrt5)^k}(a^n-b^n)^k \]
而由二项展开式可以得到：
\[(a-b)^k=C_k^0(a^n)^0(-b^n)^k+C^{k-1}_1(a^n)^1(-b^n)^{k-1}+...+C_k^k(a^n)^k(-b^n)^0 \]
那么神奇的来了：
\[\begin{equation} \begin{split} F(n)^k+F(n-1)^k&=\frac 1 {(\sqrt5)^k}((a^n-b^n)^k+(a^{n-1}-b^{n-1})^k)\\\\ &=((a^n-b^n)^k+(a^{n-1}-b^{n-1})^k)\\\\ &=C_k^0(a^n)^0(-b^n)^k+C_k^0(a^{n-1})^0(-b^{n-1})^k+...\\\\ &=C_k^0(a^0)^n((-b)^k)^n+C_k^0(a^0)^{n-1}((-b)^k)^{n-1}+...\\\\ &=(-1)^kC_k^0[(a^0b^k)^n+(a^0b^k)^{n-1}]+... \end{split} \end{equation} \]
如果继续加上\(F[n-2]+F[n-3]+...\)
会发现构成一个等比数列求和的问题。
\[(-1)^kC_k^0[(a^0b^k)^n+(a^0b^k)^{n-1}+(a^0b^k)^{n-2}+...++(a^0b^k)^{1}] \]
可以把 \(a^0b^k\) 看成公比，那么由求和公式：
\[S_n=a_1\times \frac{q^n-1}{q-1} \]
令\(a^0b^k=x\) 得：
\[\sum_{i=1}^nF(i)^k=\frac{1}{\sqrt5^k}\sum_{i=0}^k(-1)^iC_k^ix\frac{x^n-1}{x-1}(x=a^ib^{k-i}) \]
大功告成，如果看懂了这个式子，也就成功了一半了,现在我们加上\(c\) 的限制。
\[F(cn)^k=\frac 1 {(\sqrt5)^k}(a^{cn}-b^{cn})^k \]
接下来的步骤是一样的，会发现结果是：
\[\sum_{i=1}^nF(ci)^k=\frac{1}{\sqrt5^k}\sum_{i=0}^k(-1)^iC_k^ix\frac{x^n-1}{x-1}(x=(a^ib^{k-i})^c) \]
现在需要处理一下根号的问题，因为是分数取模。
首先得到\(i^2mod(1e9+9)==5\)时\(i=383008016\)
然后根据这个值得出上述\(a,b\) 的值：
#include <bits/stdc++.h>
using namespace std;
long long mod=1e9+9;
long long fastpow(long long a, long long b)
{
long long res = 1;
a%=mod;
while(b){
if(b & 1)
res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
int main(){
printf("%lld\n",(1+383008016)*fastpow(2,mod-2)%mod);
printf("%lld\n",((1-383008016)%mod+mod)%mod*fastpow(2,mod-2)%mod);
}
得到\(a=691504013,b=308495997\)
接下来就是代码的事情了，先放出来比赛时疯狂TLE的代码来帮助理解：
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e5 + 10;
const long long mod = 1e9 + 9;
long long fac[MAXN], A[MAXN], B[MAXN];
//fac数组表示阶乘，之后求C用
//A[i]表示a的i次方，B[i]同理
void Init() {
fac[0] = 1;
for (int i = 1; i < MAXN; i++)
fac[i] = fac[i - 1] * i % mod;
A[0] = B[0] = 1;
for (int i = 1; i < MAXN; i++) {
A[i] = A[i - 1] * 691504013 % mod;
//691504013表示的意义上文提到了
B[i] = B[i - 1] * 308495997 % mod;
//同上
}
}
long long fastpow(long long a, long long b) {//快速幂
long long ans = 1;
a %= mod;
while (b) {
if (b & 1)ans = ans * a % mod;
b >>= 1;
a = a * a % mod;
}
return ans;
}
long long Solve(long long n, long long c, long long k) {
long long ans = 0;
for (int i = 0; i <= k; i++) {
long long x = fastpow(A[k - i] * B[i] % mod, c);
//首先求出公比
long long C = fac[k] * fastpow(fac[k - i] * fac[i] % mod, mod - 2) % mod;
//求出组合数
long long tmp = x * (fastpow(x, n) - 1) % mod * fastpow(x - 1, mod - 2) % mod;
//求出等比数列之和
if (x == 1) tmp = n % mod;
//如果公比为1，则是一个恒等数列，特判
tmp = tmp * C % mod;
//求出最后的值
if (i & 1) ans -= tmp;
//判断是加是减
else ans += tmp;
ans %= mod;
}
long long num = fastpow(383008016, mod - 2);
//算出根号5在分母的值
ans = ans * fastpow(num, k) % mod;
//算出k次方
ans = (ans % mod + mod) % mod;
//使结果为正
return ans;
}
int main() {
int t;
long long n, k, c;
Init();
scanf("%d", &t);
while (t--) {
scanf("%lld%lld%lld", &n, &c, &k);
printf("%lld\n", Solve(n, c, k));
}
return 0;
}
然后这个代码TLE了十几次。
优化
我感觉是卡了常数，优化的话首先有一个等式：
\(a^b\%x=a^{b\%(x-1)}\%x\)
第二点，上面的代码每次都会对\(x\)取快速幂一次，这个也是很关键的一点，因为每次循坏都取。观察发现，其实只需要第一次求出
\((a^c)^0(b^c)^k\)的值，之后把这个值乘上\(\frac{a^c}{b^c}\)即可得到下一项\((a^c)^1(b^c)^{k-1}\)
那我们可以省去一个常数。
AC代码：
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e5 + 10;
const long long mod = 1e9 + 9;
long long fac[MAXN], inv[MAXN];
long long fastpow(long long a, long long b) {
long long ans = 1;
a %= mod;
while (b) {
if (b & 1)ans = ans * a % mod;
b >>= 1;
a = a * a % mod;
}
return ans;
}
void Init() {
fac[0] = inv[0] = 1;
for (int i = 1; i < MAXN; i++) {
fac[i] = fac[i - 1] * i % mod;
inv[i] = fastpow(fac[i], mod - 2);
}
}
long long Solve(long long n, long long c, long long k) {
long long ans = 0;
long long A = fastpow(691504013, c % (mod - 1)), B = fastpow(308495997, c % (mod - 1));
long long a = 1, b = fastpow(B, k), ib = fastpow(B, mod - 2);
for (int i = 0; i <= k; i++) {
long long x = a * b % mod;
long long C = fac[k] * inv[i] % mod * inv[k - i] % mod;
long long sum = x * (fastpow(x, n % (mod - 1)) - 1 + mod) % mod * fastpow(x - 1, mod - 2) % mod;
if (x == 1) sum = n % mod;
if ((k - i) & 1) ans -= sum * C % mod;
else ans += sum * C % mod;
ans %= mod;
a = a * A % mod;
b = b * ib % mod;
}
long long num = fastpow(383008016ll, mod - 2);
ans = ans * fastpow(num, k) % mod;
ans = (ans % mod + mod) % mod;
return ans;
}
int main() {
int t;
long long n, k, c;
Init();
scanf("%d", &t);
while (t--) {
scanf("%lld%lld%lld", &n, &c, &k);
printf("%lld\n", Solve(n, c, k));
}
return 0;
}