题解-[国家集训队]Crash的数字表格 / JZPTAB
前置知识：
莫比乌斯反演 </>
单组测试数据，给定 (n,m) ，求
[sumlimits_{i=1}^nsumlimits_{j=1}^moperatorname{lcm}(i,j)mod 20101009 ]
数据范围：(1le n,mle 10^7)。
作为写出了最暴力的做法的蒟蒻，来推个式子。
(nle m)，一气呵成：
[egin{split}
g(n,m)=&sumlimits_{i=1}^nsumlimits_{j=1}^moperatorname{lcm}(i,j)\
=&sumlimits_{i=1}^nsumlimits_{j=1}^mfrac{ij}{gcd(i,j)}\
=&sumlimits_{d=1}^nsumlimits_{i=1}^nsumlimits_{j=1}^mfrac{ij}{d}[gcd(i,j)=d]\
=&sumlimits_{d=1}^nsumlimits_{i=1}^{lfloorfrac nd
floor}sumlimits_{j=1}^{lfloorfrac md
floor}ijd[gcd(i,j)=1]\
=&sumlimits_{d=1}^n dsumlimits_{i=1}^{lfloorfrac nd
floor}isumlimits_{j=1}^{lfloorfrac md
floor}jsumlimits_{k|gcd(i,j)}mu(k)\
=&sumlimits_{d=1}^n dsumlimits_{k=1}^nmu(k)sumlimits_{i=1}^{lfloorfrac nd
floor}i[k|i]sumlimits_{j=1}^{lfloorfrac md
floor}j[k|j]\
=&sumlimits_{d=1}^n dsumlimits_{k=1}^nmu(k)sumlimits_{i=1}^{lfloorfrac {n}{dk}
floor}iksumlimits_{j=1}^{lfloorfrac {m}{dk}
floor}jk\
=&sumlimits_{d=1}^n dsumlimits_{k=1}^nk^2mu(k)frac{lfloorfrac{n}{dk}
floor(lfloorfrac{n}{dk}
floor+1)}{2}cdotfrac{lfloorfrac{m}{dk}
floor(lfloorfrac{m}{dk}
floor+1)}{2}\
end{split}
]
将 (x=dk) 带入：
[g(n,m)=sumlimits_{x=1}^nxcdotfrac{lfloorfrac{n}{x}
floor(lfloorfrac{n}{x}
floor+1)}{2}cdotfrac{lfloorfrac{m}{x}
floor(lfloorfrac{m}{x}
floor+1)}{2}sumlimits_{k|x}kmu(k)
]
然后筛 (mu(k)) 时顺便计算 (h(k)=kmu(k))，最后狄利克雷前缀和求 (f(k)=sumlimits_{k|x}kmu(k))。
别忘了膜拜 (20101009)，时间复杂度 (Theta(N+n))。
#include <bits/stdc++.h>
using namespace std;
//&Start
#define lng long long
#define lit long double
#define kk(i,n) "
"[i<n]
const int inf=0x3f3f3f3f;
const lng Inf=1e17;
//&Mobius
const int N=1e7;
const int mod=20101009;
bitset<N+10> np;
int mu[N+10],cnt,p[N+10],f[N+10];
void Mobius(){
f[1]=mu[1]=1;
for(int i=2;i<=N;i++){
if(!np[i]) p[++cnt]=i,mu[i]=-1;
f[i]=(mu[i]*i+mod)%mod;
for(int j=1;j<=cnt&&i*p[j]<=N;j++){
np[i*p[j]]=1;
if(i%p[j]==0){mu[i*p[j]]=0;break;}
mu[i*p[j]]=-mu[i];
}
}
for(int j=1;j<=cnt;j++)
for(int i=1;i*p[j]<=N;i++)
(f[i*p[j]]+=f[i])%=mod; //狄利克雷前缀和
}
//&Data
int n,m,ans;
int bitfun(int x){
lng res=1ll*x*f[x]%mod;
(res*=1ll*(n/x+1)*(n/x)/2%mod)%=mod;
(res*=1ll*(m/x+1)*(m/x)/2%mod)%=mod; //如上
//这个1ll不乘要爆long long，30分。
return (int)res;
}
//&Main
int main(){
Mobius();
scanf("%d%d",&n,&m);
if(n>m) swap(n,m);
for(int i=1;i<=n;i++)
(ans+=bitfun(i))%=mod;
printf("%d
",ans);
return 0;
}
祝大家学习愉快！