原題連結
考察:莫比烏斯反演
思路:
  參考的大佬的題解.這位老師總結的套路總結的很好:GO
主要為:

1. 優先提取gcd(i,j)
2. 將\(gcd(i,j) = d\)化為 \(gcd(i/d,j/d) = 1\)
3. \(\sum_{d|gcd(i,j)}^{gcd(i,j)} mob[i] = [gcd(i,j)==1]\)
4. 優先按倍數的方式列舉

解法為兩次數論分塊,時間複雜度大概在\(O(N)\) 參考題解 GO

Code

#include <iostream>
#include <cstring> 
#include <cstdio>
using namespace std;
typedef long long LL;
const int N = 1e7+10,M = 20101009;
int n,m,prime[N],cnt,mob[N],sum[N];
bool st[N];
LL p;
LL qsm(LL a,LL k,int mod)
{
	LL res=  1;
	while(k)
	{
		if(k&1) res = res*a%mod;
		a = a*a%mod;
		k>>=1;
	}
	return res;
}
void GetPrime(int n)
{
	mob[1] =1;
	for(int i=2;i<=n;i++)
	{
		if(!st[i]) prime[++cnt] = i,mob[i] = -1;
		for(int j=1;prime[j]<=n/i;j++)
		{
			st[i*prime[j]] =1;
			if(i%prime[j]==0) break;
			mob[i*prime[j]] = (-1)*mob[i];
		}
	}
	for(int i=1;i<=n;i++) sum[i] = (sum[i-1]+(LL)i*i%M*mob[i]%M)%M;
}
int calc(int x,int y)
{
	LL res = 0;
	for(int i=1,r;i<=x;i=r+1)
	{
		r = min(x/(x/i),y/(y/i));
		int a = x/i,b=y/i;
		res+=(LL)(sum[r]-sum[i-1])%M*(LL)(1+a)%M*a%M*p%M*(1+b)%M*b%M*p%M;
		res=(res+M)%M;
	}
	return res;
}
int main()
{
	scanf("%d%d",&n,&m);
	if(n>m) swap(n,m);
	GetPrime(N-1);
	LL res =0;
	p = qsm(2,M-2,M);
	for(int i=1,r;i<=n;i=r+1)
	{
		r = min(n/(n/i),m/(m/i));
		res = res+(LL)(i+r)*(r-i+1)%M*p%M*calc(n/i,m/i)%M;
		res=(res+M)%M;
	}
	printf("%lld\n",res);
	return 0;
}