Description
求
∑i=1n∑j=1mlcm(i,j)
n,m<=10^7
Solution
(注意,以下内容默认
n<=m
)
看着lcm不爽,把它变一变:
∑i=1n∑j=1mijgcd(i,j)
莫比乌斯反演常用,枚举gcd
∑d=1n1d∑i=1n∑j=1,gcd(i,j)=dmij
设
f(d)=∑i=1n∑j=1,gcd(i,j)=dmij
则
Ans=∑d=1nf(d)d
考虑反演常用,设
g(d)=∑i=1n∑j=1,d|gcd(i,j)mij
即
g(d)=∑x=1⌊nd⌋f(xd)
可以发现,
g(d)=∑i=1⌊nd⌋∑j=1⌊md⌋d2ij
g(d)=d2∑i=1⌊nd⌋∑j=1⌊md⌋ij
g(d)=d2∑i=1⌊nd⌋i∑j=1⌊md⌋j
g(d)=d2⌊nd⌋(⌊nd⌋+1)⌊md⌋(⌊nd⌋+1)4
根据反演得
f(d)=∑i=1⌊nd⌋μ(i)g(di)
f(d)=∑i=1⌊nd⌋μ(i)d2i2⌊ndi⌋(⌊ndi⌋+1)⌊mdi⌋(⌊ndi⌋+1)4
Ans=14∑d=1nd∑i=1⌊nd⌋μ(i)i2⌊ndi⌋(⌊ndi⌋+1)⌊mdi⌋(⌊ndi⌋+1)
这样发现后面的值只和 ⌊ndi⌋ 和 ⌊mdi⌋ 有关于是就可以分块了。(分块大法好!)
然后求和再用一次分块。(分块大法好!)
这样就可以 O(n) 解决问题了。(分块大法好!)
Code
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(ll i=a;i<=b;i++)
#define ll long long
#define mo 20101009
#define N 10000005
using namespace std;
int p[N],tot,mi[N],mu[N];
bool f[N];
ll n,m,ans;
ll sum(ll l,ll r) {return (r+l)*(r-l+1)/2%mo;}
ll get(ll x,ll y) {return ((x*(x+1)/2)%mo)*((y*(y+1)/2)%mo)%mo;}
int main() {
scanf("%lld%lld",&n,&m);
if (n>m) swap(n,m);mu[1]=1;
fo(i,2,n) {
if (!f[i]) p[++tot]=i,mu[i]=-1;
fo(j,1,tot) {
int k=p[j]*i;if (k>n) break;f[k]=1;
if (!(i%p[j])) break;mu[k]=-mu[i];
}
}
fo(i,1,n) mi[i]=(mi[i-1]+mu[i]*i%mo*i%mo)%mo;
ll l=1,r;
while (l<=n) {
ll nn=n / l,mm=m / l;r=min(n/nn,m / mm);
ll le=1,ri,an=0;
while (le<=nn) {
ll x=nn/le,y=mm/le;ri=min(nn/x,mm/y);
an=(an+(mi[ri]-mi[le-1])%mo*get(x,y)%mo)%mo;le=ri+1;
}
ans=(ans+sum(l,r)*an%mo)%mo;l=r+1;
}
printf("%lld",(ans+mo)%mo);
}