【莫比乌斯反演】BZOJ2154 Crash的数字表格

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

Description

求sigma lcm(x,y)，x<=n,y<=m。n,m<=1e7。

Solution

lcm没有什么直接做的好方法，用lcm=x*y/gcd转成gcd来做

就是要求sigma d*f(x/d,y/d)

f(x,y)为x和y以内gcd正好为1的对数

F为所有对数，于是有F(x,y)=x*(x+1)/2*y*(y+1)/2

f(x,y)=sigma (1<=i<=x) i*i*mu(i)*F(x/i,y/i)

f用莫比乌斯反演解决，这两个式子都套上分块优化到sqrt，于是总复杂度sqrt*sqrt=n

分块优化具体可以见上一篇blog

Code

一开始全开LL MLE了一发

然后又WA了两发，第一次是有一地方算的时候溢出

一开始为了避免MLE把prime数组/50，但其实只能/20的样子

#include<cstdio>
#include<algorithm>
#include<cstring>
#define ll long long
using namespace std;
const int maxn=1e7+,mod=; bool flag[maxn];
int prime[maxn],mu[maxn],cnt;
int sum[maxn],s[maxn];
int n,m; void getmu(){
mu[]=;
for(int i=;i<=n;i++){
if(!flag[i]){
mu[i]=-;
prime[++cnt]=i;
}
for(int j=;i*prime[j]<=n&&j<=cnt;j++){
flag[i*prime[j]]=;
if(i%prime[j]==){
mu[i*prime[j]]=;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for(int i=;i<=n;i++)
sum[i]=(sum[i-]+(ll)i*i*mu[i]%mod)%mod;
} ll F(int x,int y){
return (ll)((ll)x*(x+)/%mod)*((ll)y*(y+)/%mod)%mod;
} ll f(int x,int y){
ll ret=;
int p;
for(int i=;i<=x;i=p+){
p=min(x/(x/i),y/(y/i));
ret=(ret+(ll)(sum[p]-sum[i-])*F(x/i,y/i))%mod;
}
return ret;
} int main(){
scanf("%d%d",&n,&m);
if(n>m) swap(n,m);
for(int i=;i<=n;i++) s[i]=(s[i-]+i)%mod;
getmu(); int pos;
ll ans=;
for(int i=;i<=n;i=pos+){
pos=min(n/(n/i),m/(m/i));
ans=(ans+(ll)(s[pos]-s[i-])*f(n/i,m/i))%mod;
}
printf("%lld\n",(ans+mod)%mod);
return ;
}