題意：
給定三個數字 \(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;
}