Description

給定 \(N,C \le 10^{18}, K \le 10^5\)，對斐波那契數列 \(F\)，求 \((F_0)^K + (F_C)^K + (F_{2C})^K + ... + (F_{NC})^K\)。\(T\) 組資料。模 \(10^9+9\) 輸出。

Solution

斐波那契數列的通項為

\[F_n = d(a^n-b^n), \quad d=\frac 1 {\sqrt 5}, \quad a=\frac{1+\sqrt 5}{2}, \quad b=\frac{1-\sqrt 5}{2} \]

在模 \(10^9+9\) 下，可以暴力找到 \(\sqrt 5 = 383008016\)

於是有

\[\begin{aligned} \sum_{i=0}^N (F_{iC})^k &=d^k \sum_{i=0}^N (a^{iC}-b^{iC})^k \\ &=d^k\sum_{i=0}^N \sum_{j=0}^k(-1)^{k-i} C_k^j a^{iCj}b^{iC(k-j)} \\ & =d^k \sum_{j=0}^kC_k^j \frac {1-t^{n+1}}{1-t}, \quad t=(a^jb^{k-j})^C \end{aligned} \]

#include <bits/stdc++.h>

using namespace std;

#define int long long

const int mod = 1e9+9;
const int sq5 = 383008016;

int a,b,d,frac[100005],pwa[100005],pwb[100005],iv[100005];

int qpow(int p,signed q)
{
    return (q&1 ? p : 1) * (q ? qpow(p*p%mod, q>>1) : 1) % mod;
}

int inv(int p)
{
    if(p<=1e5) return iv[p];
    return qpow(p, mod-2);
}

void prepare()
{
    for(int i=0;i<=1e5;i++) iv[i]=qpow(i, mod-2);

    a=(1+sq5)*inv(2)%mod;
    b=(1-sq5+mod)%mod*inv(2)%mod;
    d=inv(sq5);

    for(int i=0;i<=1e5;i++) pwa[i]=qpow(a,i);
    for(int i=0;i<=1e5;i++) pwb[i]=qpow(b,i);

    frac[0]=1;
    for(int i=1;i<=1e5;i++) frac[i]=frac[i-1]*i%mod;
}

int C(int n,int m)
{
    return frac[n] * inv(frac[m]) % mod * inv(frac[n-m]) % mod;
}

int fib(int n)
{
    int ans = qpow(a,n) + qpow(b,n) + mod;
    ans %= mod;
    ans *= d;
    ans %= mod;
    return ans;
}

void work()
{
    int n,c,k;
    cin>>n>>c>>k;

    int ans=0;
    int _C=1;

    int t = pwb[k];
    t = qpow(t, c%(mod-1));

    int u = (n+1)%(mod-1);
    int tt = qpow(t,u);

    int d1 = qpow(a,c%(mod-1)), d2 = inv(qpow(b,c%(mod-1)));
    int dn1 = qpow(d1,u), dn2 = qpow(d2,u);

    for(int i=0;i<=k;i++)
    {
        if(i) t = t*d1%mod*d2%mod, tt = tt*dn1%mod*dn2%mod;

        if(i) _C = _C * (k-i+1) % mod * inv(i) % mod;

        if(t!=1) ans += ((k-i)&1 ? -1 : 1)*(_C*(tt-1+mod)%mod*inv((t-1+mod)%mod)%mod);
        else ans += ((k-i)&1 ? -1 : 1)*(_C*(n%mod+1)%mod*t%mod);
        ans %= mod;
        ans += mod;
        ans %= mod;
    }

    cout << (ans*qpow(d,k)%mod) << endl;
}

signed main()
{
    ios::sync_with_stdio(false);

    prepare();

    int t;
    cin>>t;
    while(t--)
    {
        work();
    }
}