求

\[\sum_{i=1}^N\sum_{j=1}^M(i,j)^k\]

多組詢問，但是每次的\(k\)都是不變的

先是套路

\[f(n)=\sum_{i=1}^N\sum_{j=1}^M[(i,j)=n]\]

\[F(n)=\sum_{n|d}f(d)=\left \lfloor \frac{N}{n} \right \rfloor\times\left \lfloor \frac{M}{n}\right \rfloor\]

我們要求的就是

\[Ans=\sum_{i=1}^Ni^kf(i)\]

\[=\sum_{i=1}^Ni^k\sum_{i|d}\mu(\frac{d}{i})\times \left \lfloor \frac{N}{d} \right \rfloor\times\left \lfloor \frac{M}{d}\right \rfloor\]

交換一下列舉順序

\[Ans=\sum_{d=1}\left \lfloor \frac{N}{d} \right \rfloor\times\left \lfloor \frac{M}{d}\right \rfloor\sum_{i|d}\mu(\frac{d}{i})i^k\]

我們設\(h(d)=\sum_{i|d}\mu(\frac{d}{i})i^k\)

這顯然是一個積性函式

考慮一下如何線篩這個積性函式

如果\(d\)為質數，顯然\(h(d)=d^k-1\)

發現這個形式很像尤拉函式，可以嘗試搞一下\(h(d^t)\)和\(h(d^{t-1})\)之間的關係

發現從\(d^t\)

可以直接不用考慮這個約數

於是

\[h(d^{t-1})=\sum_{i|d^{t-1}}\mu(\frac{d^{t-1}}{i})i^k\]

因為新多出來的那一個約數並不需要考慮，因為其\(\mu\)值等於0

所以\(h(d^t)\)和\(h(d^{t-1})\)的有用的約數並沒有什麼變化

所以

\[h(d^t)=\sum_{i|d^{t-1}}\mu(\frac{d^{t-1}}{i})(i\times d)^k\]

\[=d^k\sum_{i|d^{k-1}}\mu(\frac{d^{k-1}}{i})i^k=d^k\times h(d^{t-1})\]

發現這個跟尤拉函式好像啊，我們完全可以按照線性篩尤拉函式的方式來線篩\(h\)

所以線篩之後求一個字首和之後分塊就好了

#include<iostream>
#include<cstring>
#include<cstdio>
#define re register
#define maxn 5000005
#define LL long long
#define min(a,b) ((a)<(b)?(a):(b))
#define max(a,b) ((a)>(b)?(a):(b))
const LL mod=1000000007;
inline int read()
{
    char c=getchar();
    int x=0;
    while(c<'0'||c>'9') c=getchar();
    while(c>='0'&&c<='9')
        x=(x<<3)+(x<<1)+c-48,c=getchar();
    return x;
}
int p[maxn>>1],f[maxn],T,n,m;
LL pre[maxn];
int N[2005],M[2005];
inline LL quick(LL a,LL b)
{
    LL S=1;
    while(b){if(b&1) S=S*a%mod;b>>=1;a=a*a%mod;}
    return S;
}
LL K;
int main()
{
    T=read(),K=read();
    for(re int i=1;i<=T;i++) N[i]=read(),M[i]=read();
    for(re int i=1;i<=T;i++) if(N[i]>M[i]) std::swap(N[i],M[i]);
    for(re int i=1;i<=T;i++) n=max(n,M[i]);
    f[1]=1,pre[1]=1;
    for(re int i=2;i<=n;i++)
    {
        if(!f[i]) p[++p[0]]=i,pre[i]=quick(i,K)-1;
        for(re int j=1;j<=p[0]&&p[j]*i<=n;j++)
        {
            f[p[j]*i]=1;
            if(i%p[j]==0)
            {
                pre[p[j]*i]=pre[i]*(pre[p[j]]+1)%mod;
                break;
            }
            pre[p[j]*i]=pre[p[j]]*pre[i]%mod;
        }
    }
    for(re int i=1;i<=n;i++) pre[i]=(pre[i]+pre[i-1])%mod;
    for(re int t=1;t<=T;t++)
    {
        n=N[t],m=M[t];
        LL ans=0;
        for(re LL l=1,r;l<=n;l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            ans=(ans+(n/l)*(m/l)%mod*(pre[r]-pre[l-1]+mod)%mod)%mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}