[国家集训队]Crash的数字表格 / JZPTAB

题目描述

今天的数学课上，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 ;
}