题目描述
求:
∑i=1n∑j=1mlcm(i,j)
∑
i
=
1
n
∑
j
=
1
m
l
c
m
(
i
,
j
)
lcm(i,j) l c m ( i , j ) 表示i,j的最小公倍数
题解
先是一步简单的变形:
∑i=1n∑j=1mijgcd(i,j)
∑
i
=
1
n
∑
j
=
1
m
i
j
g
c
d
(
i
,
j
)
常用的套路就是枚举gcd(i,j):
∑d=1min(n,m)d∑i=1nd∑j=1md[gcd(i,j)=1](ij)
∑
d
=
1
m
i
n
(
n
,
m
)
d
∑
i
=
1
n
d
∑
j
=
1
m
d
[
g
c
d
(
i
,
j
)
=
1
]
(
i
j
)
然后其实有两种做法:
1.
∑d=1min(n,m)d∑i=1nd∑j=1md(ij)∑d′|gcd(i,j)μ(d′)
∑
d
=
1
m
i
n
(
n
,
m
)
d
∑
i
=
1
n
d
∑
j
=
1
m
d
(
i
j
)
∑
d
′
|
g
c
d
(
i
,
j
)
μ
(
d
′
)
设 F(x)表示∑xi=1i F ( x ) 表 示 ∑ i = 1 x i
∑d=1min(n,m)d∑d′=1min(nd,md)μ(d′)∗d′2∗F(ndd′)∗F(mdd′)
∑
d
=
1
m
i
n
(
n
,
m
)
d
∑
d
′
=
1
m
i
n
(
n
d
,
m
d
)
μ
(
d
′
)
∗
d
′
2
∗
F
(
n
d
d
′
)
∗
F
(
m
d
d
′
)
设dd’=T
∑T=1min(n,m)F(nT)∗F(mT)∗T∑d′|Td′μ(d′)
∑
T
=
1
m
i
n
(
n
,
m
)
F
(
n
T
)
∗
F
(
m
T
)
∗
T
∑
d
′
|
T
d
′
μ
(
d
′
)
后面那个玩意很容易筛,前面的也预处理一下,再用下求和公式就可以了。
O(n‾√) O ( n )
2.
我们直接设 f(x,y,k)=∑xi=1∑yj=1[gcd(i,j)=k]ij f ( x , y , k ) = ∑ i = 1 x ∑ j = 1 y [ g c d ( i , j ) = k ] i j
由定义可设 F(x,y,k)=∑xi=1∑yj=1[gcd(i,j)|k]ij F ( x , y , k ) = ∑ i = 1 x ∑ j = 1 y [ g c d ( i , j ) | k ] i j
(这里函数并不是说状态量多了,前两个只是范围限制,k才是真正变量)
然后原式化为:
∑d=1min(n,m)d∗f(nd,md,1)
∑
d
=
1
m
i
n
(
n
,
m
)
d
∗
f
(
n
d
,
m
d
,
1
)
∑d=1min(n,m)d∗∑1|tμ(t1)F(nd,md,t)
∑
d
=
1
m
i
n
(
n
,
m
)
d
∗
∑
1
|
t
μ
(
t
1
)
F
(
n
d
,
m
d
,
t
)
∑d=1min(n,m)d∗∑t=1min(nd,md)μ(t)F(nd,md,t)
∑
d
=
1
m
i
n
(
n
,
m
)
d
∗
∑
t
=
1
m
i
n
(
n
d
,
m
d
)
μ
(
t
)
F
(
n
d
,
m
d
,
t
)
∑d=1min(n,m)d∗∑t=1min(nd,md)μ(t)∗t2(n/d)∗(n/d+1)2∗(m/d)∗(m/d+1)2
∑
d
=
1
m
i
n
(
n
,
m
)
d
∗
∑
t
=
1
m
i
n
(
n
d
,
m
d
)
μ
(
t
)
∗
t
2
(
n
/
d
)
∗
(
n
/
d
+
1
)
2
∗
(
m
/
d
)
∗
(
m
/
d
+
1
)
2
∑d=1min(n,m)d∗(n/d)∗(n/d+1)2∗(m/d)∗(m/d+1)2∑t=1min(nd,md)μ(t)∗t2
∑
d
=
1
m
i
n
(
n
,
m
)
d
∗
(
n
/
d
)
∗
(
n
/
d
+
1
)
2
∗
(
m
/
d
)
∗
(
m
/
d
+
1
)
2
∑
t
=
1
m
i
n
(
n
d
,
m
d
)
μ
(
t
)
∗
t
2
后面那坨预处理,分块搞就可以了
复杂度
O(n‾√2)=O(n)
O
(
n
2
)
=
O
(
n
)
,但预处理常数较小,所以跑得比法1快……
法1的代码:
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
const int mod=20101009;
const int N=1e7+10;
int pri[N];bool vis[N];int mu[N];
int sum[N];int cnt=0;
int f[N];
typedef long long ll;
inline void prepare()
{
mu[1]=1;f[1]=1;
vis[1]=1;
for(register int i=2;i<N;i++)
{
if(!vis[i]) {pri[++cnt]=i;mu[i]=-1;f[i]=1-i+mod;}
for(register int j=1;j<=cnt&&(1ll*i*pri[j]<N);++j)
{
register int x=i*pri[j];
vis[x]=1;
if(i%pri[j]==0) {
mu[i]=0;f[x]=f[i];
break;
}
mu[x]=-mu[i];
f[x]=(1ll*f[i]*f[pri[j]])%mod;
}
}
sum[1]=1;
for(register int i=1;i<N;++i) {f[i]=f[i-1]+1ll*f[i]*i%mod;if(f[i]>=mod) f[i]-=mod;}// T*sigma(d*mu[d])
for(register int i=1;i<N;++i) {sum[i]=sum[i-1]+i;if(sum[i]>=mod) sum[i]-=mod;}
}
int main()
{
int n,m;prepare();
scanf("%d %d",&n,&m);register int l,r;
if(n>m) swap(n,m);register int ans=0;
for(l=1;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
register int S=f[r]-f[l-1]+mod;if(S>=mod) S-=mod;
ans+=1ll*S*sum[n/l]%mod*sum[m/l]%mod;
if(ans>=mod) ans-=mod;
}
printf("%d\n",ans);
}
法2的代码:
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#include<queue>
#include<set>
using namespace std;
typedef long long ll;
const int mod=20101009;
const int N=1e7+100;
int mu[N];ll sum[N];
int pri[N];bool vis[N];int cnt=0;
inline void prepare()
{
mu[1]=1;vis[1]=1;
for(register int i=2;i<N;i++)
{
if(!vis[i]) {pri[++cnt]=i;mu[i]=-1;}
for(register int j=1;j<=cnt&&(1ll*i*pri[j]<N);j++)
{
register int x=i*pri[j];
vis[x]=1;
if(i%pri[j]==0) {mu[x]=0;break;}
mu[x]=-mu[i];
}
}
for(register ll i=1;i<N;i++) sum[i]=(1ll*mu[i]*(i*i%mod)+sum[i-1]+mod)%mod;
}
inline ll Get(int n,int m)
{
return (((1ll*n*(n+1)/2)%mod)*((1ll*m*(m+1)/2)%mod))%mod;
}
inline ll calc(int n,int m)
{
register int l,r;
register ll res=0;
for(l=1;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
register ll S=sum[r]-sum[l-1];if(S<0) S+=mod;
res+=(1ll*S*Get(n/l,m/l))%mod;
if(res>=mod) res-=mod;
}
return res;
}
int main()
{
int n,m;prepare();
scanf("%d %d",&n,&m);
if(n>m) swap(n,m);
register ll l,r;
register ll ans=0;
for(l=1;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans+=(((l+r)*(r-l+1)/2)%mod*calc(n/l,m/l))%mod;
if(ans>=mod) ans-=mod;
}
printf("%lld\n",ans);
}