题目描述
今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时整除a和b的最小正整数。例如，LCM(6, 8) = 24。
回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张NM的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个45的表格如下：
1 2 3 4 5
2 2 6 4 10
3 6 3 12 15
4 4 12 4 20
看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod20101009的值。
输入输出格式
输入格式：
输入的第一行包含两个正整数，分别表示N和M。
输出格式：
输出一个正整数，表示表格中所有数的和mod20101009的值。
输入输出样例
输入样例#1：
4 5
输出样例#1：
122
说明
30%的数据满足N, M≤ 10^3。
70%的数据满足N, M≤ 10^5。
100%的数据满足N, M≤ 10^7。
CTM啥破玩意儿
以上只是一只SB发的牢骚==
这题是反演题
然后我不会做==
首先题目要求的是(sum_{i = 1}^{n}sum_{j=1}^{m}{frac{i * j}{gcd(i,j)}})
先让n < m
然后根据套路先考虑gcd
(sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{m{frac{i * j}{d}[gcd(i , j)== d]}})
然后继续根据常规的套路我们再把那个d提出来
(sum_{d = 1}^{n}{d}sum_{i=1}^{n/d}sum_{j=1}^{m/d}{i * j[gcd(i,j)==1]})
然后就可以设目标函数(f(x) = sum_{i=1}^{n/d}sum_{j=1}^{m/d}{i * j[gcd(i,j)==x]})
然后再设(g(x) = sum_{x|t}{f(t)})
也就是(g(x) = sum_{i=1}^{n/d}sum_{j=1}^{m/d}{i * j[x|gcd(i,j)]})
然后我们把这个x给提出来
(g(x) = x^2 sum_{i=1}^{n/dx}sum_{j=1}^{m/dx}{i * j[1|gcd(i,j)]})
然后那个([1|gcd(i,j)])没啥用了
(g(x)=x^2sum_{i=1}^{n/dx}sum_{j=1}^{m/dx}{i*j})
然后我们回到我们刚才化简的答案式子
(sum_{d = 1}^{n}{d}sum_{i=1}^{n/d}sum_{j=1}^{m/d}{i * j[gcd(i,j)==1]})
是不是可以写成(sum_{d=1}{d}f(1))
然后我们对(f(1))再反演一下
是不是就能得到(sum_{d=1}^{n}{d}sum_{i=1}^{n}{μ(i)*g(i)})
然后我们再把这个(g(i))展开
就得到了(sum_{d=1}^{n}{d}sum_{i=1}^{n}sum_{x=1}^{n/di}sum_{y=1}^{m/di}{x*y})
然后我们需要预处理出(sum_{i=1}^{n}{μ(i)*i^2})
后面的那一坨是等差数列可以O(1)计算出来
所以可以两次整除分块O(n)求解
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
# define int long long
const int M = 10000005 ;
const int mod = 20101009 ;
using namespace std ;
inline int read() {
char c = getchar() ; int x = 0 , w = 1 ;
while(c>'9'||c<'0') { if(c=='-') w = -1 ; c = getchar() ; }
while(c>='0'&&c<='9') { x = x*10+c-'0' ; c = getchar() ; }
return x*w ;
}
int mu[M] , p[M] , pnum ;
bool Notp[M] ;
int sum[M] , n , m ;
int fsum[M] , Ans ;
inline void Get_Mu() {
Notp[1] = true ; mu[1] = 1 ;
for(int i = 2 ; i <= n ; i ++) {
if(!Notp[i]) p[++pnum] = i , mu[i] = -1 ;
for(int j = 1 ; j <= pnum && i * p[j] <= n ; j ++) {
Notp[i * p[j]] = true ;
if(i % p[j] == 0) {
mu[i * p[j]] = 0 ;
break ;
}
mu[i * p[j]] = -mu[i] ;
}
}
for(int i = 1 ; i <= n ; i ++) sum[i] = (sum[i - 1] + mu[i] + mod)%mod ;
}
inline int Calc(int n , int m) {
int temp = 0 ;
for(int l = 1 , r ; l <= n ; l = r + 1) {
r = min(n / (n / l) , m / (m / l)) ;
int val = ((1 + n / l) * (n / l) / 2 % mod)%mod * ((1 + m / l) * (m / l) / 2 % mod)%mod ;
temp = (temp + (fsum[r] - fsum[l - 1] + mod)%mod * val + mod)%mod ;
}
return (temp + mod)%mod ;
}
# undef int
int main() {
# define int long long
n = read() ; m = read() ;
if(n > m) swap(n , m) ;
Get_Mu() ;
for(int i = 1 ; i <= n ; i ++) fsum[i] = (fsum[i - 1] + mu[i] * i * i%mod + mod)%mod ;
for(int l = 1 , r ; l <= n ; l = r + 1) {
r = min(n / (n / l) , m / (m / l)) ;
int val = ((l + r) * (r - l + 1LL) / 2LL + mod)%mod ;
Ans = (Ans + Calc(n / l , m / l) * val%mod)%mod ;
}
printf("%lld
",Ans) ;
return 0 ;
}