牛客练习赛84 F.牛客推荐系统开发之下班（斐波那契数列+莫比乌斯反演）

牛客练习赛84 F.牛客推荐系统开发之下班（斐波那契数列+莫比乌斯反演）

传送门：https://ac.nowcoder.com/acm/contest/11174/F
题目大意：

这是我的一个学长出的题目。。。很早之前他就让我等着这一场比赛来做他的数论题，所以这次就破例在赛中做题了。。。
题解：经典斐波那契套路题。首先你得知道斐波那契数列这样的一个性质。

$$\gcd (f_n,f_m)=f_{\gcd (n,m)}$$

有了这样一个性质，我们便可以把原式转为：

$$\sum _{i_1=1}^n\sum _{i_2=1}^n\sum _{i_3=1}^n\cdots \sum _{i_k=1}^n{f}_{\gcd (i1,i2,i3,\cdots ,i_n)}$$

接下来就是经典的莫比乌斯反演套路了。但不太一样的是，我们需要在这一步打住：

$$\sum _{d=1}^n{f}_d\sum _{t=1}^{\lfloor \frac{n}{d}\rfloor }\mu (t)\lfloor \frac{n}{dt}{\rfloor }^k$$

令 $g(n)=\sum _{t=1}^n\mu (t)\lfloor \frac{n}{t}{\rfloor }^k$ ，由于 $g(n)$ 只有 $O(\sqrt{n})$ 的不同的取值情况且 $g(n)$ 可以通过数论分块 $O(\sqrt{n})$ 求解，根据大佬的结论，暴力求解 $g(n)$ 这样的函数所有取值情况时间复杂度貌似是 $O(n^{\frac{3}{4}})$ （反正跑得飞快）。现在只需要关注 $f_d$ 的前缀和就行。还好这里还有个关于斐波那契数列前缀和 $S_n$ 的公式：

$$S_n=f_{n+2}-1$$

这下将问题转移为单点求斐波那契数列的问题了。考虑 $f_n$ 的通项公式：

$$f_n=\frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2}{)}^n-(\frac{1-\sqrt{5}}{2}{)}^n]$$

因为 $5$ 是模 $1e9+9$ 的二次剩余，说明 $\sqrt{5}$ 在模 $1e9+9$ 下存在意义。不会求二次剩余也不是很大问题，本地暴力求出即可，得到其中一个解为 $616991993$ ，即可以用 $616991993$ 替换原式中的 $\sqrt{5}$ ，那么此时原式就可以进行数论分块了。总的复杂度为 $O(n^{\frac{3}{4}})$。

#if __has_include(<bits/stdc++.h>)
#include<bits/stdc++.h>
#endif
#include<iostream>
#include<bitset>
#include<map>
#include<vector>
#include<unordered_map>
#include<algorithm>
#include<sstream>
#include<stack>
#include<queue>
#include<cmath>
using namespace std;
typedef long long ll;
typedef long double ld;
typedef pair<ll,ll> pll;
template<class T,class U>
ostream& operator<<(ostream& os,const pair<T,U>& p){
    os<<p.first<<" "<<p.second;
    return os;
}
#define debug(x) cerr<<#x<<":"<<x<<endl
#define debug2(x,y) cerr<<#x<<":"<<x<<"  "<<#y<<":"<<y<<endl
#define debug3(x,y,z) cerr<<#x<<":"<<x<<"  "<<#y<<":"<<y<<"  "<<#z<<":"<<z<<endl
#define debuga(a) for(int i=0;i<a.size();++i) debug2(i,a[i]);
#define debugar(a,b) for(int i=0;i<b;++i) debug2(i,a[i]);
struct fastio{fastio(){ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);}}fio;
#define endl '\n'
const ll N=5e6+5,M=1e9+9;
const ll C=616991993;
const ll C1=(1+C)/2,C2=(1-C)/2+M;
int mu[N],prime[N],cnt=0;
bitset<N> vis;
ll qpow(ll a,ll b){
    a%=M;
    if(a<0) a+=M;
    ll res=1;
    b%=M-1;
    if(b<0) b+=M-1;
    while(b){
        if(b&1) res=res*a%M;
        a=a*a%M;
        b>>=1;
    }
    return res;
}
void pre(){
    mu[1]=1;
    for(ll i=2;i<N;++i){
        if(!vis[i]){
            prime[cnt++]=i;
            mu[i]=-1;
        }
        for(int j=0;i*prime[j]<N;++j){
            ll t=i*prime[j];
            vis[t]=1;
            if(i%prime[j]) mu[t]=mu[i]*mu[prime[j]];
            else{
                mu[t]=0;
                break;
            }
        }
    }
    for(ll i=2;i<N;++i)  mu[i]=(mu[i-1]+mu[i])%M;
}
unordered_map<ll,ll> sm;
ll cal_sm(ll n){
    if(n<N) return mu[n];
    if(sm.count(n)) return sm[n];
    ll res=1;
    for(ll l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        res=(res-(r-l+1)*cal_sm(n/l))%M;
    }
    return sm[n]=res;
}
unordered_map<ll,ll> kmi;
void init(ll n,ll k){
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        kmi[n/l]=qpow(n/l,k);
    }
}
ll cal_g(ll n,ll k){
    ll res=0;
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        res=(res+(cal_sm(r)-cal_sm(l-1))%M*kmi[n/l]%M)%M;
    }
    return res;
}
const ll INV=qpow(C,-1);
ll cal_sf(ll n){
    if(!n) return 0;
    return ((qpow(C1,n+2)-qpow(C2,n+2))*INV%M-1)%M;
}
int main(){
    pre();
    ll n,k,res=0;
    cin>>n>>k;
    init(n,k);
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        res=(res+(cal_sf(r)-cal_sf(l-1))%M*cal_g(n/l,k)%M)%M;
    }
    if(res<0) res+=M;
    cout<<res<<endl;
}