链接：https://ac.nowcoder.com/acm/contest/549/J
来源：牛客网
来源：牛客网
题目描述
小A最近开始研究数论题了，这一次他随手写出来一个式子，∑ni=1∑mj=1gcd(i,j)2∑i=1n∑j=1mgcd(i,j)2，但是他发现他并不太会计算这个式子，你可以告诉他这个结果吗，答案可能会比较大，请模上1000000007。
输入描述:
一行两个正整数n，m一行两个正整数n，m
输出描述:
一行一个整数表示输出结果一行一个整数表示输出结果
小A最近开始研究数论题了，这一次他随手写出来一个式子，∑ni=1∑mj=1gcd(i,j)2∑i=1n∑j=1mgcd(i,j)2，但是他发现他并不太会计算这个式子，你可以告诉他这个结果吗，答案可能会比较大，请模上1000000007。
输入描述:
一行两个正整数n，m一行两个正整数n，m
输出描述:
一行一个整数表示输出结果一行一个整数表示输出结果
输入：
2 2
输出：
7
2 2
输出：
7
1=<n,m<=1e6
解题思路：这题应该算是一题莫比乌斯反演的套路题了，感觉莫比乌斯真的好难啊，学了好久感觉也没懂，这题算是它的一个简单应用。
具体可以参考博客：https://blog.sengxian.com/algorithms/mobius-inversion-formula
具体可以参考博客：https://blog.sengxian.com/algorithms/mobius-inversion-formula
莫比乌斯反演经典套路：
现在有个积性函数f(n)，设n<m，则：
于是原来的式子就变成了求f∗μ了，再用上整数分块就可以快速搞定了。
自己推演了一遍：
代码：
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#include<string>
#include<set>
#include<cmath>
#include<list>
#include<deque>
#include<cstdlib>
#include<bitset>
#include<stack>
#include<map>
#include<cstdio>
#include<queue>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
#define lson l,mid,rt<<1
#define rson mid+1,r,rt<<1|1
const int INF=0x3f3f3f3f;
const double PI=acos(-1.0);
const double eps=1e-;
const int maxn=;
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){return a/gcd(a,b)*b;}
const int mod=1e9+;
const int dir[][]={{,},{-,},{,},{,-}};
const int N=1e6+;
ll n,m,prime[N],mu[N],tot;
void getMu(){
for(int i=;i<=1e6+;i++) prime[i]=;
mu[]=;
for(int i=;i<=1e6+;i++){
if(prime[i]){
prime[++tot]=i;
mu[i]=-;
}
for(int j=;j<=tot&&prime[j]*i<=1e6+;j++){
prime[i*prime[j]]=;
if(i%prime[j]==){
mu[i*prime[j]]=;
break;
}else mu[i*prime[j]]=-mu[i];
}
}
}
int main(){
cin>>n>>m;
getMu();
ll ans=;
for(ll i=;i<=min(n,m);i++){
ll tmp=;
for(ll j=i;j<=min(n,m);j+=i){
tmp=(tmp+mu[j/i]*(n/j)*(m/j))%mod;
}
ans=(ans+tmp*i*i%mod)%mod;
}
cout<<ans<<endl;
return ;
}
整除分块优化：
#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const int maxn=1e6+;
const int mod=1e9+;
ll n,m,prime[maxn],mu[maxn],sum[maxn],tot,ans;
void getMobius(int N){
for(int i=;i<=N;i++)prime[i]=;
mu[]=;
for(int i=;i<=N;i++){
if(prime[i]){
prime[tot++]=i;
mu[i]=-;
}
for(int j=;j<tot&&i*prime[j]<=N;j++){
prime[i*prime[j]]=;
if(i%prime[j]==){
mu[i*prime[j]]=;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
}
ll solve(ll a,ll b){
ll res=;
for(int l=,r;l<=min(a,b);l=r+){
r=min(a/(a/l),b/(b/l));
res=(res+(sum[r]-sum[l-])%mod*(a/l)%mod*(b/l)%mod)%mod;
}
return res;
}
int main(){
scanf("%lld%lld",&n,&m);
if(n>m) swap(n,m);
getMobius(1e6);
sum[]=;
for(int i=;i<=1e6;i++) sum[i]=sum[i-]+mu[i];
for(ll l=,r;l<=n;l=r+){
r=min(n/(n/l),m/(m/l));
ll dd=(r*(r+)*(*r+)/-(l-)*l*(*l-)/)%mod;
ans=(ans+dd*solve(n/l,m/l)%mod)%mod;
}
printf("%lld\n",ans);
return ;
}