牛客练习赛#84 F 莫比乌斯反演+杜教筛+技巧+斐波那契数列和gcd的结论+矩阵快速幂

题意：

给出$n,k$，计算
$$\sum _{i_1=1}^n\sum _{i_2=1}^n...\sum _{i_n=1}^{n}gcd(f_{i_1},f_{i_2},...,f_{i_n})$$

Solution：

这题结论有点多，先给出来，后附第二条的证明，第一条难度太大

结论$1$：

$$gcd(f_{i_1},f_{i_2},....,f_{i_n})=f_{gcd(i_1,i_2,...,i_n)}$$

于是即计算
$$\sum _{i_1=1}^n\sum _{i_2=1}^n...\sum _{i_n=1}^n{f}_{gcd(i_1,i_2,...,i_n)}$$
显然枚举因子，这里因子为$gcd$项，即
$$\sum _{d=1}^n{f}_d\sum _{i_1=1}^n\sum _{i_2=1}^n...\sum _{i_n=1}^n[gcd(i_1,i_2,...,i_n)=d]$$
这是个莫比乌斯反演的经典形式了
$$\sum _{d=1}^n{f}_d\sum _{i_1=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{i_2=1}^{\lfloor \frac{n}{d}\rfloor }...\sum _{i_n=1}^{\lfloor \frac{n}{d}\rfloor }[gcd(i_1,i_2,...,i_n)=1]=\sum _{d=1}^n{f}_d\sum _{i_1=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{i_2=1}^{\lfloor \frac{n}{d}\rfloor }...\sum _{i_n=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{t|gcd(i_1,i_2,...,i_n)}\mu (t)$$
再次优先枚举因子，这里为$t$
$$\sum _{d=1}^n{f}_d\sum _{t=1}^{\lfloor \frac{n}{d}\rfloor }\mu (t)\lfloor \frac{n}{dt}{\rfloor }^k$$

到这里，我的思路为枚举$T=dt$，即计算
$$\sum _{T=1}^n\lfloor \frac{n}{T}{\rfloor }^k\sum _{d|T}\mu (d)f(\frac{T}{d})$$
第二个求和形式极其像莫比乌斯反演的变换形式，但很可惜，斐波那契数列不是一个积性函数，否则将是绝杀，只能作罢

不能化简的话，有一个技巧，第二个求和发现上界是一个下取整，把他令成一个函数
$$R(x)=\sum _{i=1}^x\mu (i)\lfloor \frac{x}{i}\rfloor$$
先把他代入原式我们再分析他的妙处，原式即
$$\sum _{d=1}^n{f}_{d}R(\lfloor \frac{n}{d}\rfloor )$$
首先，这个表达式有下取整，可以数论分块；其次，$R(x)$内部也有下取整，又可以数论分块，并且求$R(x)$需要$\mu$的前缀和，这恰好就是杜教筛！

解决了$R(x)$部分，那么求上式也需要$f$的前缀和，如何处理？

结论$2$：

$$\sum _{i=1}^n{f}_i=f_{n+2}-1$$

那么只需要求斐波那契的一个项即可，矩阵快速幂即可轻松做到

结论$2$证明：

利用累加法，设$S(n)$为斐波那契数列的前$n$项和：
$$\begin{gathered} f_3=f_2+f_1 \\ {f}_4=f_3+f_2 \\ ... \\ {f}_n=f_{n-1}+f_{n-2} \end{gathered}$$
对应位置全部相加，得到：
$$S(n)-f_1-f_2=(S(n)-f_1-f_n)+(S(n)-f_{n-1}-f_n)$$
$f_1=f_2=1$，合并得到
$$S(n)=2f_n+f_{n-1}-1$$
递推公式合并右式，就得到了
$$S(n)=f_{n+2}-1$$
同理是可以得到奇数项和，偶数项和的

// #include<bits/stdc++.h>
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<vector>
#include<bitset>
#include<map>
using namespace std;

using ll=long long;
const int N=2e5+5,inf=0x3fffffff;
const long long INF=0x3fffffffffffffff,mod=1e9+9;

ll qm(ll x){return (x%mod+mod)%mod;}

struct matrix
{
    ll a[2][2];
    void initial(){memset(a,0,sizeof(a));}
    void build()
    {
        initial();
        for(int i=0;i<2;i++) a[i][i]=1;
    }
    ll* operator[](int i) const {
        return (ll*)&a[i][0];
    }
    matrix operator*(const matrix& x) const
    {
        matrix ret; ret.initial();
        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++)
                for(int k=0;k<2;k++) ret[i][j]=(ret[i][j]+a[i][k]*x[k][j]%mod)%mod;
        return ret;
    }
}ans={{
    {1,1},
    {0,0}
}},tmp={{
    {1,1},
    {1,0}
}};

matrix qpow(matrix a,ll b)
{
    matrix ret,base=a; ret.build();
    while(b)
    {
        if(b&1) ret=ret*base;
        base=base*base;
        b>>=1;
    }
    return ret;
}

ll f(int x)
{
    if(x<=2)
    {
        if(x==0) return 0;
        return 1;
    }
    return qm((ans*qpow(tmp,x-2))[0][0]);
}

ll fpre(int x){return qm(f(x+2)-1);}

ll qpow(ll a,ll b)
{
    ll ret=1,base=a;
    while(b)
    {
        if(b&1) ret=ret*base%mod;
        base=base*base%mod;
        b>>=1;
    }
    return ret;
}

int mu[N],cnt,prime[N];
bool nt[N];
map<int,ll>map1;

void make_prime()
{
    mu[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!nt[i]) prime[++cnt]=i,mu[i]=mod-1;
        for(int j=1;j<=cnt&&i*prime[j]<N;j++)
        {
            nt[i*prime[j]]=true;
            if(i%prime[j]==0) break;
            else mu[i*prime[j]]=qm(1ll*mu[i]*mu[prime[j]]);
        }
    }
    for(int i=1;i<N;i++) mu[i]=qm(1ll*mu[i]+1ll*mu[i-1]);
}

ll calcmu(int x)
{
    if(x<N) return 1ll*mu[x];
    if(map1.count(x)) return map1[x];
    ll ret=1;
    for(int l=2,r;l<=x;l=r+1)
    {   
        r=x/(x/l);
        ret=qm(ret-(r-l+1)*calcmu(x/l));
    }
    return map1[x]=ret;
}

int n,k;

ll R(int x)
{
    ll ret=0;
    for(int l=1,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ret=qm(ret+qm(calcmu(r)-calcmu(l-1))*qpow(x/l,k));
    }
    return ret;
}

int main()
{
    make_prime();
    cin>>n>>k; ll res=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        res=qm(res+qm(fpre(r)-fpre(l-1))*R(n/l));
    }
    cout<<res;
    return 0;  
}