于神之怒加强版
题意即求：
\[\sum _ {i =1} ^ {n} \sum _ { j = 1 } ^ {m} \gcd (i,j) ^ k \]其中 \(k\) 是与数据组数一起给定的。
数据满足 \(T\le 2 \times 10 ^ 3\)，\(n,m,k\le 5 \times 10 ^ 6\)。
直接推导：
\[\begin{aligned} & \sum _ {i =1} ^ {n} \sum _ { j = 1 } ^ {m} \gcd (i,j) ^ k \\ = & \sum _ {d = 1} \sum _ {i =1} ^ {n} \sum _ { j = 1 } ^ {m} \left[ \gcd (i,j) = d \right] d ^ k \\ = & \sum _ {d = 1} d ^ k \sum _ {i =1} ^ {\left\lfloor n / d\right \rfloor} \sum _ { j = 1 } ^ {\left\lfloor m / d\right \rfloor} \left[ \gcd (i,j) = 1 \right] \\ = & \sum _ {d = 1} d ^ k \sum _ {i =1} ^ {\left\lfloor n / d\right \rfloor} \sum _ { j = 1 } ^ {\left\lfloor m / d\right \rfloor} \sum _ {e \mid i , e \mid j} \mu (e) \\ = & \sum _ {d = 1} d ^ k \sum _ {e=1} \mu (e) \left\lfloor \dfrac{n}{de}\right \rfloor \left\lfloor \dfrac{m}{de}\right \rfloor \\ = & \sum _ {T=1} \sum _ {d \mid T} \mu \left(\dfrac{T}{d}\right) d ^ k \left\lfloor \dfrac{n}{T}\right \rfloor \left\lfloor \dfrac{m}{T}\right \rfloor \\ =& \sum _ {T=1} \left(\sum _ {d \mid T} \mu \left(\dfrac{T}{d}\right) d ^ k \right)\left\lfloor \dfrac{n}{T}\right \rfloor \left\lfloor \dfrac{m}{T}\right \rfloor \end{aligned} \]括号内函数可以 \(O(n \ln n)\) 预处理出前缀和。
后面的是个整除分块。
时间复杂度 \(O(n\ln n + T\sqrt{n})\)。
代码：
#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
namespace Ehnaev{
inline ll read() {
ll ret=0,f=1;char ch=getchar();
while(ch<48||ch>57) {if(ch==45) f=-f;ch=getchar();}
while(ch>=48&&ch<=57) {ret=(ret<<3)+(ret<<1)+ch-48;ch=getchar();}
return ret*f;
}
inline void write(ll x) {
static char buf[22];static ll len=-1;
if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
else {do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
while(len>=0) putchar(buf[len--]);
}
}using Ehnaev::read;using Ehnaev::write;
inline void writeln(ll x) {write(x);putchar(10);}
const ll mo=1e9+7,N=5e6;
ll T,k,n,m,cnt;
ll f[N+5],sf[N+5];
ll mu[N+5],prime[N+5];
bool ff[N+5];
inline ll Qpow(ll b,ll p) {
ll r=1;while(p) {if(p&1) r=r*b%mo;b=b*b%mo;p>>=1;}return r;
}
inline void Init() {
ff[1]=1;mu[1]=1;
for(ll i=2;i<=N;i++) {
if(!ff[i]) {prime[++cnt]=i;mu[i]=mo-1;}
for(ll j=1;j<=cnt&&prime[j]*i<=N;j++) {
ff[i*prime[j]]=1;
if(i%prime[j]==0) {mu[i*prime[j]]=0;break;}
mu[i*prime[j]]=mu[i]*(mo-1)%mo;
}
}
for(ll i=1;i<=N;i++) {
ll tmp=Qpow(i,k);
for(ll j=i,c=1;j<=N;j+=i,c++) {f[j]=(f[j]+mu[c]*tmp%mo)%mo;}
}
for(ll i=1;i<=N;i++) {sf[i]=(sf[i-1]+f[i])%mo;}
}
int main() {
T=read();k=read();Init();
while(T--) {
ll ans=0;
n=read();m=read();
for(ll i=1,j;i<=n&&i<=m;i=j+1) {
j=min(n/(n/i),m/(m/i));
ans=(ans+(((sf[j]-sf[i-1]+mo)%mo)*(n/i)%mo)*(m/i)%mo)%mo;
}
writeln(ans);
}
return 0;
}