[題目鏈接] https://www.luogu.org/problemnew/show/P4449
給定n,m,k,計算
$\sum_{i=1}^n \sum_{j=1}^m \mathrm{gcd}(i,j)^k$
對1000000007取模的結果

// luogu-judger-enable-o2
/*
-----------------------
基本套路:
1.枚舉約數
2.枚舉整除分塊,觀察倍數關系
3.發現g(T)是[積性函數],且g=μ*f.此時有很好的轉移方法
  if(i為質數) g[i]=f[i]-1
  else{
     if(i為某個質數整數冪) g[p^k]=g[p^(k-1)]*f[p]+f[1]*μ[p^k]
     else g[i]=g[i/low[i]]*g[low[i]*prime[j]]   即把最小的約數都放到一起,以滿足互質
  }
4.預處理的初始化,g[1]=1.
-----------------------2019.2.15
*/
#include<bits/stdc++.h>
using namespace std;
#define int long long
typedef long long LL;
const int INF=1e9+7;
inline LL read(){
    register LL x=0,f=1;register char c=getchar();
    while(c<48||c>57){if(c==‘-‘)f=-1;c=getchar();}
    while(c>=48&&c<=57)x=(x<<3)+(x<<1)+(c&15),c=getchar();
    return f*x;
}

const int MAXN=5e6+5;
const int mod=1e9+7;

LL mu[MAXN],g[MAXN],f[MAXN],sum[MAXN],prime[MAXN],low[MAXN];
bool vis[MAXN];
int T,n,m,k;

inline int qpow(int a,int b){
    LL res=1;
    while(b){
        if(b&1) (res*=a)%=mod;
        (a*=a)%=mod;
        b>>=1;
    }
    return res;
}

inline void init(int n){
    mu[1]=1;
    g[1]=1;//
    for(int i=1;i<=n;i++)
        f[i]=qpow(i,k);
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            prime[++prime[0]]=i;
            mu[i]=-1;
            low[i]=i;
            g[i]=(f[i]-1)%mod;//i為質數的轉移
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;j++){
            vis[i*prime[j]]=true;
            if(i%prime[j]==0){
                low[i*prime[j]]=low[i]*prime[j];
                if(low[i]==i)  //整次冪情況
                    g[i*prime[j]]=(g[i]*f[prime[j]])%mod;// +f[1]*mu[i*prime[j]]
                else
                    g[i*prime[j]]=(g[i/low[i]]*g[low[i]*prime[j]])%mod;
                break;
            }
            else{
                g[i*prime[j]]=(g[i]*g[prime[j]])%mod;
                low[i*prime[j]]=prime[j];//每個數只會由它最小的約數更新一次
                mu[i*prime[j]]=-mu[i];
            }
        }
    }
    for(int i=1;i<=n;i++)
        sum[i]=(sum[i-1]+g[i])%mod;
}

signed main(){
    //freopen("4449.in","r",stdin);
    T=read(),k=read();
    init(5e6);
    while(T--){
        n=read(),m=read();
        if(n>m) swap(n,m);
        LL ans=0;
        for(int l=1,r;l<=n;l=r+1){
            r=min(n/(n/l),m/(m/l));
            (ans+=(n/l)*(m/l)%mod*(sum[r]-sum[l-1]+mod)%mod)%=mod;
        }
        printf("%lld\n",ans);
    }
}
 

P4449 於神之怒加強版 （莫比烏斯反演）