bzoj3561 DZY Loves Math VI

​​http://www.elijahqi.win/archives/3456​​​
Description

给定正整数n,m。求

Input

一行两个整数n,m。
Output

一个整数，为答案模1000000007后的值。
Sample Input

5 4
Sample Output

424
HINT

数据规模：

1<=n,m<=500000，共有3组数据。

Source

By Jcvb
∑ni=1∑mj=1lcm(i,j)gcd(i,j)=∑nd=1∑⌊nd⌋i=1∑⌊md⌋j=1(ijd)d⋅[gcd(i,j)=1]=∑nd=1∑⌊nd⌋i=1∑⌊md⌋j=1(ijd)d⋅∑p|i,p|jμ(p)=∑nd=1dd∑⌊nd⌋p=1μ(p)∑⌊npd⌋i=1∑⌊mpd⌋j=1(ijp2)d=∑nd=1dd∑⌊nd⌋p=1μ(p)p2d∑⌊npd⌋i=1id∑⌊mpd⌋j=1jd ∑ i = 1 n ∑ j = 1 m l c m ( i , j ) gcd ( i , j ) = ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ ( i j d ) d ⋅ [ gcd ( i , j ) = 1 ] = ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ ( i j d ) d ⋅ ∑ p | i , p | j μ ( p ) = ∑ d = 1 n d d ∑ p = 1 ⌊ n d ⌋ μ ( p ) ∑ i = 1 ⌊ n p d ⌋ ∑ j = 1 ⌊ m p d ⌋ ( i j p 2 ) d = ∑ d = 1 n d d ∑ p = 1 ⌊ n d ⌋ μ ( p ) p 2 d ∑ i = 1 ⌊ n p d ⌋ i d ∑ j = 1 ⌊ m p d ⌋ j d

#include<cstdio>
#include<algorithm>
#define ll long long
using namespace std;
const int mod=1e9+7;
inline void inc(int &x,int v){x=x+v>=mod?x+v-mod:x+v;}
inline void dec(int &x,int v){x=x-v<0?x-v+mod:x-v;}
inline int ksm(ll b,int t){static ll tmp;
    for (tmp=1;t;b=b*b%mod,t>>=1) if (t&1) tmp=tmp*b%mod;return tmp; 
}
const int N=5e5+10;
bool not_prime[N];
int mu[N],prime[N],tot,n,m,sum[N],jc[N];
inline void init(){
    mu[1]=1;
    for (int i=2;i<=5e5;++i){
        if (!not_prime[i]) prime[++tot]=i,mu[i]=-1;
        for (int j=1;prime[j]*i<=5e5;++j){
            not_prime[prime[j]*i]=1;
            if (i%prime[j]==0){
                mu[prime[j]*i]=0;break;
            }else mu[prime[j]*i]=-mu[i];
        }   
    }
}
int main(){
    freopen("bzoj3561.in","r",stdin);
    scanf("%d%d",&n,&m);if (n>m) swap(n,m);init();int ans=0;
    for (int i=1;i<=m;++i) jc[i]=1;
    for (int d=1,tmp;d<=n;++d){
        for (int i=1;i*d<=m;++i) jc[i]=(ll)jc[i]*i%mod,sum[i]=(sum[i-1]+jc[i])%mod;tmp=0;
        for (int i=1;i*d<=n;++i){
            if (mu[i]==1) inc(tmp,(ll)ksm(i,d<<1)*sum[n/d/i]%mod*sum[m/d/i]%mod);
            else if (mu[i]==-1) dec(tmp,(ll)ksm(i,d<<1)*sum[n/d/i]%mod*sum[m/d/i]%mod);
        }inc(ans,(ll)tmp*ksm(d,d)%mod);
    }printf("%d\n",ans);
    return 0;
}