题目描述
今天的数学课上, 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^2\sum_{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\n",Ans) ;
- return 0 ;
- }
来源: http://www.bubuko.com/infodetail-2752487.html