【bzoj】P4407于神之怒加强版（莫比乌斯反演）

　　题目链接

　　套路一般的枚举$gcd(i,j)=w$。设$min(n,m)=top$,则有

　　$\sum\limits_{i=1}^{n}\sum\limits_{j=1}{m}gcd(i,j)$

　　$=\sum\limits_{w=1}^{top}w^{k}\sum\limits_{w|i}^{n}\sum\limits_{w|j,(i,j)=w}^{m}1$

　　我们设$f(w)=\sum\limits_{w|i}^{n}\sum\limits_{w|j,(i,j)=w}^{m}1$

　　$F(w)=\sum\limits_{w|i}^{n}\sum\limits_{w|j}^{m}1$

　　则有$F(w)=\sum\limits_{w|d}f(d)$

　　然后就根据莫比乌斯反演公式

　　$f(w)=\sum\limits_{w|d}\mu(\frac{d}{w})F(w)$

　　然后容易想到$F(w)=\frac{n}{w}\frac{m}{w}$

　　然后就有了原式

　　$=\sum\limits_{w=1}^{top}w^{k}\sum\limits_{d=1}^{\frac{top}{w}}\mu(d)\frac{n}{wd}\frac{m}{wd}$

　　然后……Timelimitexceed，我就看题解了。

　　枚举t=wd。然后数论分块乱搞。

　　讲道理我如果抗住题解的诱惑自己推的话还是可以推出来的qwq。

　　

#include<cstdio>
#include<cctype>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#define maxn 6000010
#define mod 1000000007
using namespace std;
inline long long read(){
    long long num=0,f=1;
    char ch=getchar();
    while(!isdigit(ch)){
        if(ch=='-')    f=-1;
        ch=getchar();
    }
    while(isdigit(ch)){
        num=num*10+ch-'0';
        ch=getchar();
    }
    return num*f;
}

int miu[maxn];
bool s[maxn];
int prime[1000010],tot;
long long mul[maxn];
long long sum[maxn];

long long Pow(long long n,int i,long long p){
    long long ret=1;
    while(i){
        if(i&1)    ret=(ret*n)%p;
        n=(n*n)%p;
        i>>=1;
    }
    return ret;
}

int main(){
    int T=read(),L=read();
    miu[1]=mul[1]=sum[1]=1;
    for(int i=2;i<maxn;++i){
        if(!s[i]){
            prime[++tot]=i;
            mul[tot]=Pow(i,L,mod);
            sum[i]=mul[tot]-1;
        }
        for(int j=1;j<=tot&&1LL*i*prime[j]<=maxn;++j){
            s[i*prime[j]]=1;
            if(i%prime[j])    sum[i*prime[j]]=(1ll*sum[i]*sum[prime[j]])%mod;
            else{
                sum[i*prime[j]]=(1ll*sum[i]*mul[j])%mod;
                break;
            }
        }
    }
    for(int i=2;i<maxn;++i){
        sum[i]+=sum[i-1];
        sum[i]%=mod;
    }
    while(T--){
        int n=read(),m=read();
        int top=min(n,m);
        int x=1; long long ans=0;
        if(n>m)    swap(n,m);
        while(x<=top){
            int y=min(n/(n/x),m/(m/x));
            ans+=1ll*(n/x)*(m/x)%mod*(sum[y]-sum[x-1])%mod;
            ans%=mod;
            x=y+1;
        }
        printf("%lld\n",(ans+mod)%mod);
    }
    return 0;
}

/*
1 2
3 3
*/