[HDU6755]Fibonacci Sum

最新推荐文章于 2024-08-26 15:30:00 发布

原创最新推荐文章于 2024-08-26 15:30:00 发布 · 1.6k 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #数学

数学同时被 3 个专栏收录

3 篇文章

订阅专栏

二项式定理

1 篇文章

订阅专栏

分块幂

1 篇文章

订阅专栏

本文深入解析了斐波那契数列的高级求和算法，利用矩阵快速幂和模意义下的快速幂优化，实现了在大规模数据范围内的高效计算。文章详细介绍了如何通过预处理和分块幂技巧，将复杂度降低到O(∑klog⁡P)，并提供了完整的代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题意

定义Fibonacci 数列：

$F_0=0,F_1=1\\ F_n=F_{n-1}+F_{n-2}\ (n>1)$

给定整数 $N, c, k$ ，求
$\sum_{n=0}^NF_{nc}^k\mod{10^9+9}$
$N,c\le10^{18},k\le10^5$ .

题解

根据Fibonacci数列通项公式 $\displaystyle F_n=\frac{1}{\sqrt5}\Big[(\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n\Big]$ ，有
$\begin{aligned} \sum_{n=0}^NF_{nc}^k&=\sum_{n=0}^N\Big(\frac{1}{\sqrt5}\Big)^k\Big[(\frac{1+\sqrt5}{2})^{nc}-(\frac{1-\sqrt5}{2})^{nc}\Big]^k,\quad a_c=(\frac{1+\sqrt5}{2})^c,b_c=(\frac{1-\sqrt5}{2})^c\\ &=\Big(\frac{1}{\sqrt5}\Big)^k\sum_{n=0}^N\sum_{i=0}^k(-1)^i{k\choose i}(a_c^{k-i}b_c^i)^n,\quad t_i=a_c^{k-i}b_c^i\\ &=\Big(\frac{1}{\sqrt5}\Big)^k\sum_{i=0}^k(-1)^i{k\choose i}\sum_{n=0}^Nt_i^n\\ &=\Big(\frac{1}{\sqrt5}\Big)^k\sum_{i=0}^k(-1)^i{k\choose i}\cdot\left\{\begin{aligned} \frac{t_i^{N+1}-1}{t_i-1},&\quad t_i\neq 1\\ N+1,&\quad t_i=1 \end{aligned}\right. \end{aligned}$

暴力预处理出 $383008016^2\equiv616991993^2\equiv5\pmod{10^9+9}$ ，任取一个作为 $\sqrt5$ 模意义下的值即可。
再求出 $\displaystyle a\equiv\frac{1+\sqrt5}{2}\pmod{10^9+9},b\equiv\frac{1-\sqrt5}{2}\pmod{10^9+9)}$ 。
计算快速幂的时候用 $a^b\equiv a^{b\mod{P-1}}\pmod{P}$ 降幂优化。
$t_i$ 可以通过 $\displaystyle t_i=t_{i-1}\cdot(\frac ba)^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);
}

预处理分块幂 $\displaystyle a^n=a^{S\lfloor\frac{n}{S}\rfloor+(n\mod{S})}=a^{S\lfloor\frac{n}{S}\rfloor}\cdot a^{n\mod{S}}=f_{a,1}(\lfloor\frac{n}{S}\rfloor)\cdot f_{a,2}(n\mod{S})$ ，取 $S=2^{15},S^2>10^9+9$ ，将取模和整除转化为位运算。这样快速幂就可以 $O (1)$ 求出。
假设 $t_i\neq 1$ ，则有
$\sum_{n=0}^NF_{nc}^k=\Big(\frac{1}{\sqrt5}\Big)^k\frac{1}{\prod\limits_{i=0}^k(t_i-1)}\sum_{i=0}^k(-1)^i{k\choose i}(t_i^{N+1}-1)\prod_{j\neq i}(t_i-1)$

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

这样替换就可以在求 $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;
}

~~卡常害人~~