BZOJ_P3561 DZY Loves Math VI(数论+莫比乌斯反演)

BZOJ传送门
Time Limit: 10 Sec Memory Limit: 256 MB
Submit: 90 Solved: 66
[Submit][Status][Discuss]
Description
给定正整数n,m。求

Input
一行两个整数n,m。

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

Sample Input
5 4

Sample Output
424

HINT
数据规模：
1<=n,m<=500000，共有3组数据。

Source
By Jcvb

#include<cstdio>
#include<iostream>
using namespace std;
#define N 500005
#define P 1000000007
typedef long long ll;
bool b[N];int mu[N],p[N];
int n,m,cnt,ans;int sum[N],a[N];
inline int Pow(int a,int b,int res=1){while(b){if(b&1) res=(ll)res*a%P;a=(ll)a*a%P,b>>=1;}return res;}
void GetMu(){
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!b[i]) mu[i]=-1,p[cnt++]=i;
        for(int j=0;j<cnt&&i*p[j]<=n;j++){
            b[i*p[j]]=1;
            if(i%p[j]) mu[i*p[j]]=-mu[i];
            else break;
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);if(n<m) swap(n,m);GetMu();
    for(int i=1;i<=n;i++) a[i]=1;
    for(int d=1;d<=m;d++){
        int x=Pow(d,d),y=0,limm=m/d,limn=n/d;
        for(int j=1;j<=limn;j++) a[j]=(ll)a[j]*j%P,sum[j]=(sum[j-1]+a[j])%P;
        for(int k=1;k<=limm;k++)if(mu[k])
            y=((ll)a[k]*a[k]%P*sum[limm/k]%P*sum[limn/k]%P*mu[k]+y+P)%P;
        ans=(ans+(ll)x*y%P)%P;
    }
    printf("%d",ans);
    return 0;
}