BZOJ 3561 DZY Loves Math VI（莫比乌斯反演）

Description

给定正整数n,m，求∑i=1n∑j=1mlcm(i,j)gcd(i,j)

Input

一行两个整数n,m，1≤n,m≤5e5

Output

一个整数，为答案模109+7后的值

Sample Input

5 4

Sample Output

424

Solution

不妨假设n≤m

∑i=1n∑j=1mlcm(i,j)gcd(i,j)====∑i=1n∑j=1m(igcd(i,j)jgcd(i,j)gcd(i,j))gcd(i,j)∑d=1ndd∑i=1nd∑j=1mdidjd[(i,j)=1]∑d=1ndd∑i=1nd∑j=1mdidjd∑k=1ndμ(k)∑d=1ndd∑k=1ndμ(k)k2d∑i=1ndkid∑j=1mdkjd

枚举d，预处理Sum(i)=∑i=1ndid，总时间复杂度O(nlogn)

Code

#include<cstdio>#include<iostream>#include<cstring>#include<algorithm>#include<cmath>#include<vector>#include<queue>#include<map>#include<set>#include<ctime>using namespace std;typedef long long ll;typedef pair<int,int>P;const int INF=0x3f3f3f3f,maxn=500001;#define mod 1000000007int mu[maxn],p[maxn],mark[maxn],res=0;void init(int n=5e5){    mu[1]=1;    for(int i=2;i<=n;i++)    {        if(!mark[i])p[res++]=i,mu[i]=-1;        for(int j=0;j<res&&i*p[j]<=n;j++)        {            mark[i*p[j]]=1;            if(i%p[j])mu[i*p[j]]=-mu[i];            else            {                mu[i*p[j]]=0;                break;            }        }    }}int Pow(int a,int b){    int ans=1;    while(b)    {        if(b&1)ans=(ll)ans*a%mod;        a=(ll)a*a%mod;        b>>=1;    }    return ans;}void add(int &x,int y){    x=x+y>=mod?x+y-mod:x+y;}int n,m,a[maxn],s[maxn];int main(){    init();    scanf("%d%d",&n,&m);    if(n>m)swap(n,m);    int ans=0;    for(int i=1;i<=m;i++)a[i]=1;    for(int d=1;d<=n;d++)    {        for(int i=1;i*d<=m;i++)a[i]=(ll)a[i]*i%mod,s[i]=(s[i-1]+a[i])%mod;        int temp=0;        for(int i=1;i*d<=n;i++)        {            if(mu[i]==1)add(temp,(ll)Pow(i,2*d)*s[n/d/i]%mod*s[m/d/i]%mod);            else if(mu[i]==-1)add(temp,mod-(ll)Pow(i,2*d)*s[n/d/i]%mod*s[m/d/i]%mod);        }        add(ans,(ll)Pow(d,d)*temp%mod);     }    printf("%d\n",ans);    return 0;}