2020 Multi-University Training Contest 1 - E. Fibonacci Sum

HDU 6755 - Fibonacci Sum

題意

定義\(F_i\)為斐波那契數列第\(i\)項

給定三個正整數\(N,C,K\)

求\((F_0)^k+(F_C)^k+(F_{2C})^k+...+(F_{NC})^k\)的值，結果對\(10^9+9\)取模

資料範圍

\(1≤T≤200\)

\(1≤N,C≤10^{18} , 1≤K≤10^5\)

解

《Fibonacci數列的冪和》ACdreamers這篇部落格的想法幫助很大，可以說是這道題的變種

可能是原題的題：ZOJ3774 - Power of Fibonacci

可以不需要知道的東西是：\(5\)是模數\(10^9+9\)的二次剩餘

斐波那契數列的通項公式為\(F_n=\frac{1}{\sqrt 5}[(\frac{1+\sqrt 5}{2})^n-(\frac{1-\sqrt 5}{2})^n]\)

根據逆元和二次剩餘的轉化，可以使用\(x^2≡5\ mod\ 10^9+9\)中的\(x\)來代替\(\sqrt 5\)

故可直接計算得到下列常數

\(\frac{1}{\sqrt5}≡276601605\ ,\ \frac{1+\sqrt 5}{2}≡691504013\ ,\ \frac{1-\sqrt 5}{2}≡308495997\ (mod\ 10^9+9)\)

設\(D=\frac{1}{\sqrt 5}\ ,\ a=\frac{1+\sqrt 5}{2}\ ,\ b=\frac{1-\sqrt 5}{2}\)

則原式可看作\(F_n=D(a^n-b^n)\)

\(F_n^k=D^k(a^n-b^n)^k\)

發現\(D^k\)是個常數，不會隨項數變化而變化，故可以放在最後乘上

現在我們只觀察\((a^n-b^n)^k\)

根據二項式定理，將該式子二項式展開可得

\((a^n-b^n)^k=C_k^0(a^n)^k(-b^n)^0+C_k^1(a^n)^{k-1}(-b^n)^1+C_k^2(a^n)^{k-2}(-b^n)^2+...+C_k^r(a^n)^{k-r}(-b^n)^r+...+C_k^k(a^n)^0(-b^n)^k\)

將負數項提出，得到

\((a^n-b^n)^k=C_k^0(a^n)^k(b^n)^0(-1)^0+C_k^1(a^n)^{k-1}(b^n)^1(-1)^1+...+C_k^r(a^n)^{k-r}(b^n)^r(-1)^r+...+C_k^k(a^n)^0(b^n)^k(-1)^k\)

故最終的式子為

\(F_n^k=D^k[C_k^0(a^n)^k(b^n)^0(-1)^0+C_k^1(a^n)^{k-1}(b^n)^1(-1)^1+...+C_k^r(a^n)^{k-r}(b^n)^r(-1)^r+...+C_k^k(a^n)^0(b^n)^k(-1)^k]\)

根據題意中需要我們求的式子，發現下標每項遞增一個常數\(C\)

類似的，我們可以得到

\(F_{2n}^k=D^k[C_k^0(a^{2n})^k(b^{2n})^0(-1)^0+C_k^1(a^{2n})^{k-1}(b^{2n})^1(-1)^1+...+C_k^r(a^{2n})^{k-r}(b^{2n})^r(-1)^r+...+C_k^k(a^{2n})^0(b^{2n})^k(-1)^k]\)

\(F_{3n}^k=D^k[C_k^0(a^{3n})^k(b^{3n})^0(-1)^0+C_k^1(a^{3n})^{k-1}(b^{3n})^1(-1)^1+...+C_k^r(a^{3n})^{k-r}(b^{3n})^r(-1)^r+...+C_k^k(a^{3n})^0(b^{3n})^k(-1)^k]\)

容易發現，當\(k\)固定時，每一項的\(C_k^r\)和\((-1)^r\)是相同的，而不同的項為

\[(a^n)^{k-r}(b^n)^r\\ (a^{2n})^{k-r}(b^{2n})^r\\ (a^{3n})^{k-r}(b^{3n})^r \]

將指數化開可得

\[(a^n)^{k-r}(b^n)^r\\ (a^n)^{(k-r)+(k-r)}(b^n)^{r+r}\\ (a^n)^{(k-r)+(k-r)+(k-r)}(b^n)^{r+r+r} \]

很容易能夠看出來，這正是一個等比數列

首項為\((a^n)^{k-r}(b^n)^r\)，公比也為\((a^n)^{k-r}(b^n)^r\)

故可以通過等比數列求和公式\(S=\frac{a_1(1-q^n)}{1-q}\)快速求出

有了這個性質，本題便能列舉\(k\)後直接計算答案

優化一

同時也可以發現，設相鄰兩項為

\[C_k^r(a^{n})^{k-r}(b^{n})^r(-1)^r\\ C_k^{r+1}(a^{n})^{k-(r+1)}(b^{n})^{r+1}(-1)^{r+1} \]

後一項可以由前一項乘上\(-\frac{b^n}{a^n}\)得來

故該部分改成遞推可以再減少幾次快速冪的時間

優化二

發現本題資料達到\(10^{18}\)，且模數為質數，故可以稍微改動下快速冪部分，將尤拉降冪應用進去

簡單來說，尤拉降冪即\(a^n≡a^{n\%(mod-1)+(mod-1)}\)，但存在著使用限制

優化三

對於組合數的求法，注意到題目中\(k\)的範圍僅有\(10^5\)，故可以\(O(n)\)預處理出\(10^5\)內的階乘\(fac[i]\)以及階乘逆元\(inv[i]\)

藉助組合數公式\(C_m^n=\frac{m!}{(m-n)!n!}\)直接獲得組合數

遍歷複雜度為\(O(k)\)，快速冪複雜度為\(O(logn)\)，總時間複雜度\(O(Tklogn)\)

完整程式

還是得卡常數，注意優化計算過程

最後跑出來是(1716ms/2000ms)，實在不會優化了

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

const int maxn=100005;
const ll mod=1e9+9;

ll fac[maxn+10],inv[maxn+10];

ll qpow(ll a,ll n)
{
    ll r=1;
    a%=mod; //防止自乘越界
    if(n>mod)
        n=n%(mod-1)+mod-1; //尤拉降冪
    while(n)
    {
        if(n&1)
            r=(r*a)%mod;
        n>>=1;
        a=(a*a)%mod;
    }
    return r;
}

void init()
{
    fac[0]=1;
    for(int i=1;i<=maxn;i++)
        fac[i]=fac[i-1]*i%mod; //求出階乘
    inv[maxn]=qpow(fac[maxn],mod-2); //求出最大項階乘的逆元
    for(int i=maxn-1;i>=0;i--)
        inv[i]=inv[i+1]*(i+1)%mod; //遞推得到所有階乘的逆元
}

inline ll getC(ll m,ll n) //組合數，下m上n
{
    return fac[m]*inv[m-n]%mod*inv[n]%mod; //據組合數公式直接計算組合數
}

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    init(); //初始化階乘與階乘逆元
    int T;cin>>T;
    while(T--)
    {
        ll n,c,k,ans=0;
        cin>>n>>c>>k;
        ll tmp1=qpow(691504013,c),tmp2=qpow(308495997,c);
        ll r1=qpow(tmp1,k),r2=1; //an與bn的初始值
        ll rinv=qpow(tmp1,mod-2); //每次an項需要除以a，預處理逆元
        for(int r=0;r<=k;r++)
        {
            ll t=r1*r2%mod,tmp; //兩項相乘
            if(t==1) //如果等於1，說明等比數列首項與公比均為1，直接計算
                tmp=n%mod;
            else //否則套上等比數列求和公式
                tmp=t*(qpow(t,n)-1)%mod*qpow(t-1,mod-2)%mod;
            tmp=tmp*getC(k,r)%mod; //乘上組合數
            if(r&1) //判斷該項的正負性
                ans-=tmp;
            else
                ans+=tmp;
            r1=r1*rinv%mod; //an項每次除以a
            r2=r2*tmp2%mod; //bn項每次乘上b
        }
        ans=(ans%mod+mod)%mod*qpow(276601605,k)%mod; //注意將ans計算至[0,mod)的範圍內先，最後乘上D^k
        cout<<ans<<'\n';
    }
    return 0;
}