HDU 多校训练赛第一场1005 Fibonacci Sum

Posted 2020-12-18 jianjinxiani

题意：

? 给定三个数字 (n,c,k) ，求以下式子

? (sum_{i=0}^nF(ic)^k\%(10^9+9))

? 其中(F(x))为斐波那契数列第(x)项。

? (1leq n,cleq10^{18},1leq 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)

? 那么神奇的来了：

? (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}]+...)

? 如果继续加上(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 imes 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)^kxfrac{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)^kxfrac{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
",(1+383008016)*fastpow(2,mod-2)%mod);
    printf("%lld
",((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
", 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
", Solve(n, c, k));
    }
    return 0;
}