2020 Multi-University Training Contest 1 1005 ( HDU 6755 ) Fibonacci Sum

\(Fibonacci\)数列的通项公式：

\[F_n=\frac{1}{\sqrt5}\bigg[(\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n\bigg] \]

令：

\[a=\frac{1+\sqrt5}{2} \]

\[b=\frac{1-\sqrt5}{2} \]

则：

\[F_n=\frac{1}{\sqrt5}\big(a^n-b^n\big) \]

\[F_n^k=\bigg(\frac{1}{\sqrt5}\bigg)^k\big(a^n-b^n\big)^k \]

\[F_n^k=\bigg(\frac{1}{\sqrt5}\bigg)^k\displaystyle\sum_{r=0}^{k}\big(-1\big)^rC_k^r\big(a^{k-r}b^{r}\big)^n \]

根据题意：

\[F_c^k=\bigg(\frac{1}{\sqrt5}\bigg)^k\displaystyle\sum_{r=0}^{k}\big(-1\big)^rC_k^r\big(a^{k-r}b^{r}\big)^c \]

\[F_{2c}^k=\bigg(\frac{1}{\sqrt5}\bigg)^k\displaystyle\sum_{r=0}^{k}\big(-1\big)^rC_k^r\big(a^{k-r}b^{r}\big)^{2c} \]

\[... \]

\[F_{nc}^k=\bigg(\frac{1}{\sqrt5}\bigg)^k\displaystyle\sum_{r=0}^{k}\big(-1\big)^rC_k^r\big(a^{k-r}b^{r}\big)^{nc} \]

令：

\[t_r=(a^{k-r}b^r)^c \]

则：

\[F_c^k+F_{2c}^k+...+F_{nc}^k=\bigg(\frac{1}{\sqrt5}\bigg)^k\displaystyle\sum_{r=0}^{k}\big(-1\big)^rC_k^r\frac{t_r(t_r^n-1)}{t_r-1} \]

最后处理根式，在本题中\(5\)是模\(1e9+9\)的二次剩余，所以有\(x^2\equiv5(mod\) \(1000000009)\)，解得\(x=383008016\)或\(616991993\)，任取一解代替\(\sqrt5\)
若取\(x=383008016\)，则有：

\[a=\frac{1+\sqrt5}{2}\equiv(1+\sqrt5)*inv(2)\equiv(1+383008016)*500000005\equiv691504013(mod1000000009) \]

\[b=\frac{1-\sqrt5}{2}\equiv(1-\sqrt5)*inv(2)\equiv(1-383008016)*500000005\equiv308495997(mod1000000009) \]

特判\(t_r=1\)的情况

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD = 1000000009;
const int INVS5 = 276601605; // 1/sqrt(5)
const int A = 691504013;
const int B = 308495997;
const int maxn = 100010;
int base, last;

int qpow(int x, int y)
{
    int ans = 1;
    while (y)
    {
        if (y & 1)
            ans = (ll)x * ans % MOD;
        x = (ll)x * x % MOD;
        y >>= 1;
    }
    return ans;
}

int fac[maxn], invf[maxn], a[maxn], b[maxn];

void init(int n)
{
    fac[0] = a[0] = b[0] = 1;

    for (int i = 1; i <= n; i++)
    {
        fac[i] = (ll)fac[i - 1] * i % MOD;
        a[i] = (ll)a[i - 1] * A % MOD;
        b[i] = (ll)b[i - 1] * B % MOD;
    }
    invf[n] = qpow(fac[n], MOD - 2);
    for (int i = n - 1; i >= 0; i--)
        invf[i] = (ll)invf[i + 1] * (i + 1) % MOD;
}

inline int C(int n, int m){ // n >= m >= 0
    return n < m || m < 0 ? 0 : (ll)fac[n] * invf[m] % MOD * invf[n - m] % MOD;
}
int solve(ll n, ll c, int k)
{
    int ans = 0;
    for (int i = 0; i <= k; i++)
    {
        int tr=last,q;
        last = (ll)last * base % MOD;
        if (tr == 1)
            q = n % MOD;
        else
            q = (ll)tr * (qpow(tr, n % (MOD - 1)) - 1) % MOD * qpow(tr - 1, MOD - 2) % MOD;
        int res = (ll)C(k, i) * q % MOD;
        if (i & 1)
            ans -= res;
        else
            ans += res;
        ans %= MOD;
    }
    ans = (ll)ans * qpow(INVS5, k) % MOD;
    ans = (ans + MOD) % MOD;
    return ans;
}

int t, k;
ll n, c;
int main()
{
    init(100000);
    scanf("%d", &t);
    while (t--)
    {
        scanf("%lld%lld%d", &n, &c, &k);
        base = (ll)qpow(B, c % (MOD - 1)) * qpow(qpow(A, c % (MOD - 1)), MOD - 2) % MOD;
        last = qpow(a[k], c % (MOD - 1));
        printf("%d\n", solve(n, c, k));
    }
    return 0;
}

\(refer to:\) https://blog.csdn.net/acdreamers/article/details/23039571