• 2020杭电多校 1E / HDU 6755


    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<<'
    ';
        }
        return 0;
    }
    
  • 相关阅读:
    pgsql查询优化之模糊查询
    两款不同应用场景的Wpf分页控件
    C# 客户端程序调用外部程序的三种实现
    WPF DataGrid自动生成序号
    WPF控件自适应屏幕
    NOPI实现导入导出泛型List,支持自定义列
    WPF蒙板弹窗
    定时发送邮件
    数据库基础知识介绍(MySQL)
    动态数据透视表
  • 原文地址:https://www.cnblogs.com/stelayuri/p/13357775.html
Copyright © 2020-2023  润新知