ZOJ-3774 Power of Fibonacci——等比数列求和&&等价替换

题目

求 $displaystyle sum_{i=1}^n F_i^k$，（$1 leq nleq 10^{18},1 leq  kleq 10^5$），答案对 $10^9+9$ 取模。

分析

 将通项公式 $fib_i = frac{1}{sqrt{5}} ((frac{1 + sqrt{5}}{2})^i - (frac{1 - sqrt{5}}{2})^i)$ 代入，可以得到

$$egin{align*} S & = (frac{1}{sqrt{5}})^k sumlimits_{i=1}^n ((frac{1 + sqrt{5}}{2})^i - (frac{1 - sqrt{5}}{2}) ^ i)^k \ & = (frac{1}{sqrt{5}})^k sumlimits_{i=1}^n sumlimits_{j=0}^k (-1)^{k-j} inom{k}{j}(frac{1 + sqrt{5}}{2})^{ij} (frac{1 - sqrt{5}}{2})^{i(k-j)} \ &= (frac{1}{sqrt{5}})^k sumlimits_{j=0}^k (-1)^{k-j} inom{k}{j} sumlimits_{i=1}^n [(frac{1 + sqrt{5}}{2})^{j} (frac{1 - sqrt{5}}{2})^{k-j}]^i \ &= (frac{1}{sqrt{5}})^k sumlimits_{j=0}^k (-1)^{k-j} inom{k}{j} (frac{t^{n+1} - t}{t-1})
end{align*}$$

因为 $x^2 equiv 5(mod p)$，最终结果不含 $sqrt 5$， 肯定是被平方了，所以可以用 $x$ 代替 $sqrt 5$。

因为5在模 $10^9+9$意义下有二次剩余，所以 $sqrt 5$ 有实际意义，那么我们可以从小到大枚举 $j$，后面那一部分是等比数列求和，注意特判公比为1.

如果5在某些模数下没有二次剩余，因为 $a sqrt 5 + b% 在上述需要的运算（加、减、乘、除和幂）中是封闭的，所有我们可以用 $pair(a, b)$ 表示 $a sqrt 5 + b$，并进行运算。

#include<bits/stdc++.h>
using namespace std;

#define int long long
inline int read(){
    int a = 0;
    char c = getchar();
    bool f = 0;
    while(!isdigit(c) && c != EOF){
        if(c == '-')
            f = 1;
        c = getchar();
    }
    if(c == EOF)
        exit(0);
    while(isdigit(c)){
        a = a * 10 + c - 48;
        c = getchar();
    }
    return f ? -a : a;
}

const int MOD = 1e9 + 9 , INV2 = (MOD + 1) >> 1;  //2*(p+1)/2=1
int n, k;

template < class T >
T poww(T a , int b){
    T times = 1;
    while(b){
        if(b & 1) times = times * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return times;
}

struct PII{
    int st , nd;
    PII(int _st = 0 , int _nd = 0) : st(_st) , nd(_nd){}
    PII operator =(int b){return *this = PII(b , 0);}
    bool operator !=(PII a){return st != a.st || nd != a.nd;}
};
PII operator +(PII a , PII b){return PII((a.st + b.st) % MOD , (a.nd + b.nd) % MOD);}
PII operator -(PII a , PII b){return PII((a.st + MOD - b.st) % MOD , (a.nd + MOD - b.nd) % MOD);}
PII operator *(PII a , PII b){return PII((a.st * b.st + 5 * a.nd * b.nd) % MOD , (a.st * b.nd + a.nd * b.st) % MOD);}
PII operator *(PII a , int b){return PII(a.st * b % MOD , a.nd * b % MOD);}
PII operator %(PII a , int b){return a;}
PII operator /(PII a , PII b){return a * PII(b.st , MOD - b.nd) * poww((b.st * b.st - 5 * b.nd * b.nd % MOD + MOD) % MOD , MOD - 2);}

int solve(int x , int k){
    PII all(0 , 0);
    int C = 1 , sgn = poww(MOD - 1 , k);
    for(int j = 0 ; j <= k ; ++j){
        PII cur = poww(PII(INV2 , INV2) , j) * poww(PII(INV2 , MOD - INV2) , k - j);
        if(cur != PII(1 , 0))
            all = all + (poww(cur , x + 1) - cur) / (cur - PII(1 , 0)) * sgn * C;
        else
            all = all + PII(x % MOD , 0) * sgn * C;
        C = C * (k - j) % MOD * poww(j + 1 , MOD - 2) % MOD;
        sgn = sgn * (MOD - 1) % MOD;
    }
    all = all * poww(PII(0 , poww(5LL , MOD - 2)) , k);  //模板要求poww的参数类型相同
    return all.st;
}

signed main(){
    for(int T = read() ; T ; --T){
        n = read(); k = read();
        printf("%lld
" , solve(n , k));
    }
    return 0;
}

参考链接：https://www.cnblogs.com/Itst/p/10735935.html