题目
求：
\[\large \sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k
\]
\(1\le n,m,k\le 5\times 10^6,1\le t\le 2000\)
分析
稍微推一下吧：
\[\large
\begin{split}
&\ \ \ \ \sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k\\
&=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\vert \gcd(i,j)} d^k[\gcd(\frac id,\frac jd)=1]\\
&=\sum_{d=1}^{\min(n,m)}d^k\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[\gcd(i,j)=1]\\
&=\sum_{d=1}^{\min(n,m)}d^k\sum_{k=1}^{\lfloor\frac{\min(n,m)}{d} \rfloor}\mu(k)\lfloor\frac n{kd}\rfloor\lfloor\frac m{kd}\rfloor
\end{split}
\]
后面那个东西其实就是一个数论分块，前面的就是 \(id^k\) ，根据狄利克雷卷积的性质，显然可以也是个积性函数，再加上求 \(d\) 为质数的复杂度在 \(\ln\) 以内，所以可以线性筛。
代码
#include<bits/stdc++.h>
using namespace std;
template <typename T>
inline void read(T &x){
x=0;bool f=false;char ch=getchar();
while(!isdigit(ch)){f|=ch=='-';ch=getchar();}
while(isdigit(ch)){x=x*10+(ch^48);ch=getchar();}
x=f?-x:x;
return ;
}
template <typename T>
inline void write(T x){
if(x<0) x=-x,putchar('-');
if(x>9) write(x/10);
putchar(x%10^48);
return ;
}
#define ll long long
#define ull unsigned long long
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define dep(i,y,x) for(int i=(y);i>=(x);i--)
const int N=5e6+5,M=2e5+5,MOD=1e9+7;
int n,m,k,prime[N],cnt;
int f[N],Pow[N],pre[N],low[N],Ans;
bool vis[N];
inline int inc(int x,int y){x+=y;return x>=MOD?x-MOD:x;}
inline int dec(int x,int y){x-=y;return x<0?x+MOD:x;}
inline void incc(int &x,int y){x+=y;if(x>=MOD) x-=MOD;}
inline void decc(int &x,int y){x-=y;if(x<0) x+=MOD;}
inline int QuickPow(int x,int y){
int res=1;
while(y){
if(y&1) res=1ll*res*x%MOD;
x=1ll*x*x%MOD;
y>>=1;
}
return res;
}
inline void GetPrimes(int d){
vis[1]=true;f[1]=low[1]=pre[1]=1;//对1进行定义
for(ll i=2;i<=d;i++){
if(!vis[i]) prime[++cnt]=low[i]=i,Pow[cnt]=QuickPow(i,k),f[i]=dec(Pow[cnt],1);//对质数进行定义
pre[i]=inc(pre[i-1],f[i]);
for(ll j=1;j<=cnt&&i*prime[j]<=d;j++){
vis[i*prime[j]]=true;
if(i%prime[j]==0){
low[i*prime[j]]=low[i]*prime[j];
if (low[i]==i) f[i*prime[j]]=1ll*f[i]*Pow[j]%MOD;//对质数的若干次幂进行定义(一般由f[i]递推);
else f[i*prime[j]]=1ll*f[i/low[i]]*f[low[i]*prime[j]]%MOD;//把当前质因子的所有项剃掉，然后把新的加回来
break;
}
low[i*prime[j]]=prime[j];
f[i*prime[j]]=1ll*f[i]*f[prime[j]]%MOD;//通过 f(n),f(m) 计算 f(nm)
}
}
return ;
}
signed main(){
int t;read(t),read(k);GetPrimes(5e6);
while(t--){
read(n),read(m);
if(n>m) swap(n,m);Ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
incc(Ans,1ll*(n/l)*(m/l)%MOD*(pre[r]-pre[l-1]+MOD)%MOD);
}
write(Ans);putchar('\n');
}
return 0;
}