[HDU6755]Fibonacci Sum

题意

定义Fibonacci 数列：

$$\begin{gathered} F_0=0,F_1=1 \\ {F}_n=F_{n-1}+F_{n-2}(n>1) \end{gathered}$$

给定整数$N,c,k$，求
$$\sum _{n=0}^N{F}_{nc}^k\mathrm{m}\mathrm{o}\mathrm{d}10^9+9$$
$N,c\le 10^{18},k\le 10^5$.

---

题解

根据Fibonacci数列通项公式$F_n=\frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2}{)}^n-(\frac{1-\sqrt{5}}{2}{)}^n]$，有
$$\begin{array}{rl}\sum _{n=0}^N{F}_{nc}^k& =\sum _{n=0}^N\left(\frac{1}{\sqrt{5}}{)}^k\right[(\frac{1+\sqrt{5}}{2}{)}^{nc}-(\frac{1-\sqrt{5}}{2}{)}^{nc}{]}^k,a_c=(\frac{1+\sqrt{5}}{2}{)}^c,b_c=(\frac{1-\sqrt{5}}{2}{)}^c\\ & =(\frac{1}{\sqrt{5}}{)}^k\sum _{n=0}^N\sum _{i=0}^k(-1{)}^i\left(\genfrac{}{}{0ex}{}{k}{i}\right)(a_c^{k-i}{b}_c^i{)}^n,t_i=a_c^{k-i}{b}_c^i\\ & =(\frac{1}{\sqrt{5}}{)}^k\sum _{i=0}^k(-1{)}^i\left(\genfrac{}{}{0ex}{}{k}{i}\right)\sum _{n=0}^N{t}_i^n\\ & =(\frac{1}{\sqrt{5}}{)}^k\sum _{i=0}^k(-1{)}^i\left(\genfrac{}{}{0ex}{}{k}{i}\right)\cdot \{\begin{array}{rl}\frac{t_i^{N+1}-1}{t_i-1},& t_i\ne 1\\ N+1,& t_i=1\end{array}\end{array}$$

暴力预处理出$383008016^2\equiv 616991993^2\equiv 5(\mathrm{m}\mathrm{o}\mathrm{d}10^9+9)$ ，任取一个作为$\sqrt{5}$模意义下的值即可。
再求出$a\equiv \frac{1+\sqrt{5}}{2}(\mathrm{m}\mathrm{o}\mathrm{d}10^9+9),b\equiv \frac{1-\sqrt{5}}{2}(\mathrm{m}\mathrm{o}\mathrm{d}10^9+9))$。
计算快速幂的时候用$a^b\equiv a^{b\mathrm{m}\mathrm{o}\mathrm{d}P-1}(\mathrm{m}\mathrm{o}\mathrm{d}P)$降幂优化。
$t_i$可以通过$t_i=t_{i-1}\cdot (\frac{b}{a}{)}^c$递推计算。
时间复杂度$O(\sum k\log P)$

const int P = 1e9 + 9, sqrt5 = 383008016, invsqrt5 = 276601605, A = 691504013,
          B = 308495997;
inline int FastExp(int a, ll b) {
    int x = 1;
    for (b %= P - 1; b; b >>= 1, a = (ll)a * a % P)
        if (b & 1)
            x = (ll)x * a % P;
    return x;
}
inline void Solve() {
    ll Ans = 0;
    int t = FastExp(FastExp(A, k), c),
        q = FastExp((ll)B * FastExp(A, P - 2) % P, c);
    for (int i = 0; i <= k; ++i, t = (ll)t * q % P)
        if (t == 1)
            Ans += ((i & 1) ? -1 : 1) * (ll)C(k, i) * ((n + 1) % P) % P;
        else
            Ans += ((i & 1) ? -1 : 1) * (ll)C(k, i) * (FastExp(t, n + 1) - 1) % P * FastExp(t - 1, P - 2) % P;

    Ans = (Ans % P + P) % P * FastExp(invsqrt5, k) % P;
    printf("%lld\n", Ans);
}

预处理分块幂$a^n=a^{S\lfloor \frac{n}{S}\rfloor +(n\mathrm{m}\mathrm{o}\mathrm{d}S)}=a^{S\lfloor \frac{n}{S}\rfloor }\cdot a^{n\mathrm{m}\mathrm{o}\mathrm{d}S}=f_{a,1}(\lfloor \frac{n}{S}\rfloor )\cdot f_{a,2}(n\mathrm{m}\mathrm{o}\mathrm{d}S)$，取$S=2^{15},S^2>10^9+9$，将取模和整除转化为位运算。这样快速幂就可以$O(1)$求出。
假设$t_i\ne 1$，则有
$$\sum _{n=0}^N{F}_{nc}^k=(\frac{1}{\sqrt{5}}{)}^k\frac{1}{\prod _{i=0}^k(t_i-1)}\sum _{i=0}^k(-1{)}^i\left(\genfrac{}{}{0ex}{}{k}{i}\right)(t_i^{N+1}-1)\prod _{j\ne i}(t_i-1)$$

预处理出$t_i-1$的前缀积和后缀积，就可以只需要算一次逆元，再特判一下存在$t_i=1$的情况即可。
注意到$ab=-1$，于是有
$$a^{k-i}{b}^i=\{\begin{array}{rl}(-1{)}^i{a}^{k-2i},& i\le \frac{k}{2}\\ (-1{)}^{k-i}{b}^{2i-k},& i>\frac{k}{2}\end{array}$$

这样替换就可以在求$t_i$和$t_i^{N+1}$时少做几次乘法和取模运算。
时间复杂度$O(\sqrt{P}+\sum k)$

#include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 5, S = (1 << 15) - 1, P = 1e9 + 9, sqrt5 = 383008016,
          invsqrt5 = 276601605, A = 691504013, B = 308495997;
typedef long long ll;
typedef int arr[N];
/*-----------------------------------------------------------------*/
int k;
arr g, pre, suf, fac = {1}, ifac, f_1[2] = {{1}, {1}}, f_2[2] = {{1}, {1}};
ll n, c;
inline int C(int n, int m) {
    return n >= m ? (ll)fac[n] * ifac[n - m] % P * ifac[m] % P : 0;
}
inline int FastExp(int a, int b) {
    int x = 1;
    for (; b; b >>= 1, a = (ll)a * a % P)
        if (b & 1)
            x = (ll)x * a % P;
    return x;
}
inline int BlockExp(int t, int b) {//t=0表示底数为A;t=1表示底数为B
    return (ll)f_1[t][b >> 15] * f_2[t][b & S] % P;
}
inline ll Sign(ll x) { return x & 1 ? -1 : 1; }
inline void Solve() {
    //prepare
    int Sign1 = c & 1, Sign2 = (c & 1) * ((n + 1) & 1);//记录几个特殊指数的奇偶性
    c %= P - 1;
    ll Ans = 0, n1c = (n + 1) % (P - 1) * c % (P - 1);
    for (int i = 0; i <= (k - 1) / 2; ++i)
        g[i] = Sign(Sign1 * i & 1) * BlockExp(0, c * (k - 2 * i) % (P - 1)) - 1;
    if (!(k & 1))
        g[k / 2] = Sign((Sign1 * k / 2) & 1) == 1 ? 1 : -2;
    for (int i = k / 2 + 1; i <= k; ++i)
        g[i] = Sign(Sign1 * (k - i) & 1) * BlockExp(1, c * (2 * i - k) % (P - 1)) - 1;
    pre[0] = 1, suf[k] = 1;
    for (int i = 1; i <= k; ++i)
        pre[i] = (ll)pre[i - 1] * g[i - 1] % P;
    for (int i = k - 1; i >= 0; --i)
        suf[i] = (ll)suf[i + 1] * g[i + 1] % P;
    //calculate answer
    for (int i = 0; i <= (k - 1) / 2; ++i)
        Ans += Sign(i) * C(k, i) * (Sign(Sign2 * i) * BlockExp(0, n1c * (k - 2 * i) % (P - 1)) - 1) %  P * pre[i] % P * suf[i] % P;
    if (!(k & 1)) {
        int t = Sign(Sign1 * k / 2);
        Ans += Sign(k / 2) * C(k, k / 2) * (t == 1 ? (n + 1) % P : (n + 1) & 1) % P * pre[k] % P * g[k] % P;
    }
    for (int i = k / 2 + 1; i <= k; ++i)
        Ans += Sign(i) * C(k, i) * (Sign(Sign2 * (k - i)) * BlockExp(1, n1c * (2 * i - k) % (P - 1)) - 1) % P * pre[i] % P * suf[i] % P;
    Ans = (Ans % P + P) % P * FastExp(invsqrt5, k) % P * FastExp((ll)pre[k] * g[k] % P, P - 2) % P;
    printf("%lld\n", (Ans % P + P) % P);
}
int main() {
    for (int i = 1; i <= 1e5; ++i)
        fac[i] = (ll)fac[i - 1] * i % P;
    ifac[100000] = FastExp(fac[100000], P - 2);
    for (int i = 1e5; i >= 1; --i)
        ifac[i - 1] = (ll)ifac[i] * i % P;
    for (int i = 1; i <= S; ++i)
        f_2[0][i] = (ll)f_2[0][i - 1] * A % P,
        f_2[1][i] = (ll)f_2[1][i - 1] * B % P;
    f_1[0][1] = (ll)f_2[0][S] * A % P;
    f_1[1][1] = (ll)f_2[1][S] * B % P;
    for (int t = 0; t <= 1; ++t)
        for (int i = 2; i <= S; ++i)
            f_1[t][i] = (ll)f_1[t][i - 1] * f_1[t][1] % P;
    scanf("%*d");
    while (~scanf("%lld%lld%d", &n, &c, &k))
        Solve();
    return 0;
}

卡常害人