Given (n,c,k), compute (\sum_{i=0}^n f_{ic}^k) modulo (p=10^9+9), where (f_i) is the (i)-th Fibonacci number. Constraitns: (n,c\in[1,10^{18}]), (k\le 10^5), (T\le 200), time limit 3 seconds.

Use the closed-form expression for Fibonacci numbers: (f_i = \frac{1}{\sqrt5}\left(\alpha^i - \beta^i\right)), where (\alpha = \frac{1+\sqrt5}{2}), (\beta = \frac{1-\sqrt5}{2}). In modulo (p), (\sqrt5 \equiv 383008016).

Substitute into the sum: [ ans = \sum_{i=0}^n \left(\frac{1}{\sqrt5}(\alpha^{ic} - \beta^{ic})\right)^k = \left(\frac{1}{\sqrt5}\right)^k \sum_{i=0}^n (\alpha^{ic} - \beta^{ic})^k ]

Expand using the binomial theorem: [ ans = \left(\frac{1}{\sqrt5}\right)^k \sum_{i=0}^n \sum_{j=0}^k \binom{k}{j} \alpha^{icj} (-\beta^{ic})^{k-j} = \left(\frac{1}{\sqrt5}\right)^k \sum_{i=0}^n \sum_{j=0}^k (-1)^{k-j} \binom{k}{j} \alpha^{icj} \beta^{ic(k-j)} ]

Swap the order of summation: [ ans = \left(\frac{1}{\sqrt5}\right)^k \sum_{j=0}^k (-1)^{k-j} \binom{k}{j} \sum_{i=0}^n (\alpha^{jc} \beta^{(k-j)c})^i ]

The inner sum is a geometric series with ratio (r_j = \alpha^{jc} \beta^{(k-j)c}). If (r_j = 1), the sum equals (n+1); otherwise it equals (\frac{r_j^{n+1}-1}{r_j-1}).

Naively, each term requires several modular expoenntiations, leading to (O(Tk \log p)) with a high constant. Optimize by precomputing recurrence relations for the ratio and the last term. As (j) increases, (r_j) and (r_j^{n+1}) satisfy: [ r_{j+1} = r_j \cdot \frac{\alpha^c}{\beta^c}, \quad r_{j+1}^{n+1} = r_j^{n+1} \cdot \left(\frac{\alpha^c}{\beta^c}\right)^{n+1} ] Thus, compute (r_0 = \beta^{kc}) and (r_0^{n+1} = \beta^{kc(n+1)}), then multiply by (\frac{\alpha^c}{\beta^c}) and its ((n+1))-th power respectively for each subsequent (j). This reduces exponentiations to (O(k)) per test case.

Further optimizaton: compute the inverses of all ((r_j-1)) linearly. Let (d_j = r_j - 1) for (j=0,\dots,k). Compute the product (P = \prod_{j=0}^k d_j) and its inverse (P^{-1}). Then each (d_j^{-1} = P^{-1} \cdot \prod_{i<j} d_i \cdot \prod_{i>j} d_i). Precompute prefix and suffix products to get each inverse in (O(1)) after initial (O(k)) work. This avoids (O(k \log p)) inverse computations.

Implementation details:

• Preprocess factorials up to (k) for binomial coefficients.
• Compute (\alpha, \beta, \sqrt5^{-1}) modulo (p).
• For each test case, compute (\texttt{base} = \beta^{kc}), (\texttt{base_end} = \beta^{kc(n+1)}), (\texttt{ratio} = \alpha^c \cdot \beta^{-c}), (\texttt{ratio_end} = (\alpha^c \cdot \beta^{-c})^{n+1}).
• Loop (j=0) to (k), updating (r = \texttt{base}) and (r^{n+1} = \texttt{base_end}) multiplicatively.
• Compute (d = r-1), use linear inverse if needed.
• Term contribution: (\binom{k}{j} (-1)^{k-j} \cdot \frac{r^{n+1}-1}{r-1}) or (\binom{k}{j} (-1)^{k-j} (n+1)) when (r=1).
• Multiply final result by ((1/\sqrt5)^k).

Complexity: (O(k + \log p)) per test case for the optimized version, with (O(\log p)) for initial exponentiations.

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mod=1000000009,sqrt5=383008016;
int qpow(int x,int y){
    int res=1;
    while(y){
        if(y&1) res=res*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return res;
}
int inv(int x){return qpow(x,mod-2);}
const int K=100000;
int a,b;
int n,c,k;
int fact[K+1],factinv[K+1];
int comb(int x,int y){return fact[x]*factinv[y]%mod*factinv[x-y]%mod;}
void mian(){
    cin>>n>>c>>k;
    int ans=0;
    int base=qpow(qpow(b,k),c);
    int base_end=qpow(qpow(qpow(b,k),c),n);
    int mul=qpow(a,c)*inv(qpow(b,c))%mod;
    int mul_end=qpow(qpow(a,c),n)*inv(qpow(qpow(b,c),n))%mod;
    int cur=base, cur_end=base_end;
    int d_arr[K+2];
    for(int i=0;i<=k;i++){
        d_arr[i]=( (cur-1)%mod + mod ) %mod;
        if(d_arr[i]==0) d_arr[i]=1;
        cur=cur*mul%mod;
    }
    cur=base;
    // prefix and suffix products for linear inverse
    int pref[K+2]={1}, suff[K+3]={1};
    for(int i=1;i<=k+1;i++) pref[i]=pref[i-1]*d_arr[i-1]%mod;
    for(int i=k;i>=0;i--) suff[i]=suff[i+1]*d_arr[i]%mod;
    int prod_inv=inv(pref[k+1]);
    for(int i=0;i<=k;i++){
        int d = (cur-1)%mod;
        d = (d+mod)%mod;
        int inv_d;
        if(d==0){
            inv_d=0;
        }else{
            inv_d = prod_inv * pref[i] %mod * suff[i+1] %mod;
        }
        int term;
        if(cur==1){
            term = (n+1)%mod;
        }else{
            term = (cur_end*cur%mod - 1 + mod)%mod * inv_d %mod;
        }
        int sign = ((k-i)%2==0 ? 1 : -1);
        ans = (ans + comb(k,i) * sign * term) %mod;
        cur = cur * mul %mod;
        cur_end = cur_end * mul_end %mod;
    }
    int inv_sqrt5 = inv(sqrt5);
    ans = ans * qpow(inv_sqrt5, k) %mod;
    ans = (ans+mod)%mod;
    cout<<ans<<"\n";
}
signed main(){
    a=(1+sqrt5)*inv(2)%mod;
    b=(1-sqrt5+mod)*inv(2)%mod;
    fact[0]=factinv[0]=1;
    for(int i=1;i<=K;i++){
        fact[i]=fact[i-1]*i%mod;
        factinv[i]=inv(fact[i]);
    }
    ios::sync_with_stdio(false);
    cin.tie(0);
    int testnum;
    cin>>testnum;
    while(testnum--) mian();
    return 0;
}