【BZOJ 4407】于神之怒加强版（莫比乌斯反演+线性筛）

传送门

简单莫反后可以得到

$ans=sum_{T=1}^nfrac n Tfrac m Tf(T)$
其中$f=mu*Id_k$
然后线性筛$f$即可
对于$f(p^t)=-p^{(t-1)*K}+p^{tK}=f(p^{t-1})*p^{K}$

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define re register
#define pb push_back
#define pii pair<int,int>
#define ll long long
#define fi first
#define se second
#define bg begin
cs int RLEN=1<<20|1;
inline char gc(){
    static char ibuf[RLEN],*ib,*ob;
    (ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
    return (ib==ob)?EOF:*ib++;
}
inline int read(){
    char ch=gc();
    int res=0;bool f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
    return f?res:-res;
}
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
cs int mod=1e9+7;
inline int add(int a,int b){return (a+=b)>=mod?(a-mod):a;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return (r>=mod)?(r%mod):r;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline void Dec(int &a,int b){a-=b,a+=a>>31&mod;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=(r>=mod)?(r%mod):r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
inline int fix(int x){return (x<0)?x+mod:x;}
cs int N=5000005;
int f[N],pr[N],pw[N],tot;
bitset<N>vis;
int n,m,K;
inline void init(cs int len=N-5){
    f[1]=1;
    for(int i=2;i<=len;i++){
        if(!vis[i])pr[++tot]=i,pw[i]=ksm(i,K),f[i]=dec(pw[i],1);
        for(int j=1;j<=tot&&i*pr[j]<=len;j++){
            int p=i*pr[j];
            vis[p]=1;
            if(i%pr[j]==0){
                f[p]=mul(f[i],pw[pr[j]]);break;
            }
            f[p]=mul(f[i],f[pr[j]]);
        }
    }
    for(int i=1;i<=len;i++)Add(f[i],f[i-1]);
}
inline void solve(){
    n=read(),m=read();
    if(n<m)swap(n,m);
    int ret=0;
    for(int i=1,j;i<=m;i=j+1){
        j=min(n/(n/i),m/(m/i));
        Add(ret,mul(n/i,mul(m/i,dec(f[j],f[i-1]))));
    }
    cout<<ret<<'
';
}
int main(){
    #ifdef Stargazer
    freopen("lx.in","r",stdin);
    #endif
    int T=read();K=read();
    init();
    while(T--)solve();  
}