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 1000000007
int 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;
}