正题

题目链接:https://ac.nowcoder.com/acm/contest/11174/F

---

题目大意

给出$n,k$求
$$\sum _{i_1=1}^n\sum _{i_2=1}^n...\sum _{i_k=1}^{n}gcd(f_{i_1},f_{i_2},...,f_{i_k})$$
对$10^9+9$取模

其中$f_i$表示斐波那契数列的第$i$项

$1\le n,k\le 10^8$

---

解题思路

看上去就很莫比乌斯反演，首先把$f$提出来然后直接上莫反就是
$$\sum _{i=1}^n{f}_i\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }\mu (j)\lfloor \frac{n}{ij}{\rfloor }^k$$
这个式子其实就可以直接做了，外面$f_i$是一个整除分块，然后里面的式子也是一个整除分块。需要的算法是一个矩阵乘法求斐波那契的前缀和（当然因为$\sqrt{5}$在$10^9+9$下有二次剩余可以直接用特制方程求）还有一个杜教筛求$\mu$的前缀和。

因为整除分块套整除分块的复杂度是$O(n^{\frac{3}{4}})$，然后最里面套个杜教筛复杂度不会大很多所以能过。

然后有一个大佬的做法是把上面那个式子化一下
$$\sum _{i=1}^n\lfloor \frac{n}{i}{\rfloor }^k\sum _{d|i}{f}_d\mu (\frac{i}{d})$$
（相当于外面枚举$i\times j$然后枚举它的约数）

然后因为$\mu \ast I=ϵ$，所以$f\ast \mu \ast I=f\ast ϵ=f$。所以如果在能快速求$f$的前缀和的情况下求$f\ast \mu$的前缀和可以用杜教筛来做。

然后上面那个式子整除分块一下就好了

代码的写法是第一种

---

code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#define ll long long
using namespace std;
const ll N=5e6+10,S=3,P=1e9+9;
struct Matrix{ll a[S][S];
}c,f,am;
ll n,k,ans,cnt,mu[N],pri[N/10];
bool v[N];map<ll ,ll >mp;
Matrix operator*(Matrix &a,Matrix &b){memset(c.a,0,sizeof(c.a));for(ll i=0;i<S;i++)for(ll j=0;j<S;j++)for(ll k=0;k<S;k++)(c.a[i][j]+=a.a[i][k]*b.a[k][j]%P)%=P;return c;
}
ll power(ll x,ll b){ll ans=1;while(b){if(b&1)ans=ans*x%P;x=x*x%P;b>>=1;}return ans;
}
ll Fbi(ll b){memset(f.a,0,sizeof(f.a));am.a[0][0]=am.a[0][2]=0;am.a[0][1]=1;f.a[2][2]=1;f.a[0][1]=1;f.a[1][1]=1;f.a[1][0]=1;f.a[1][2]=1;while(b){if(b&1)am=am*f;f=f*f;b>>=1; }return am.a[0][2];
}
void Prime(){mu[1]=1;for(ll i=2;i<N;i++){if(!v[i])pri[++cnt]=i,mu[i]=-1;for(ll j=1;j<=cnt&&i*pri[j]<N;j++){v[i*pri[j]]=1;if(i%pri[j]==0)break;mu[i*pri[j]]=-mu[i];}}for(ll i=1;i<N;i++)mu[i]+=mu[i-1];return;
}
ll GetMu(ll n){if(n<N)return mu[n];if(mp[n])return mp[n];ll ans=1;for(ll l=2,r;l<=n;l=r+1)r=n/(n/l),ans=ans-(r-l+1)*GetMu(n/l);return mp[n]=ans;
}
ll g(ll n){ll ans=0;for(ll l=1,r;l<=n;l=r+1){r=n/(n/l);ll tmp=GetMu(r)-GetMu(l-1);ans=(ans+tmp*power(n/l,k)%P)%P;}return ans;
}
signed main()
{Prime();scanf("%lld%lld",&n,&k);for(ll l=1,r;l<=n;l=r+1){r=n/(n/l);ll tmp=Fbi(r)-Fbi(l-1);ans=(ans+tmp*g(n/l)%P)%P;}printf("%lld\n",(ans+P)%P);return 0;
}