Crash的数字表格
Time Limit: 20 Sec Memory Limit: 259 MB[Submit][Status][Discuss]
Description
今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下： 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只想知道表格里所有数的和mod 20101009的值。
Input
输入的第一行包含两个正整数，分别表示N和M。
Output
输出一个正整数，表示表格中所有数的和mod 20101009的值。
Sample Input
4 5
Sample Output
122
HINT
100%的数据满足N, M ≤ 10^7。
Main idea
求
Solution
先推一波式子：
然后我们只要求出了 f 就可以解决这个问题了。
然后我们就可以莫比乌斯反演出解了。
Code
1 #include<iostream>
2 #include<string>
3 #include<algorithm>
4 #include<cstdio>
5 #include<cstring>
6 #include<cstdlib>
7 #include<cmath>
8 using namespace std;
9 typedef long long s64;
10
11 const int ONE = 10000005;
12 const int MOD = 20101009;
13
14 int T;
15 int n,m;
16 bool isp[ONE];
17 int prime[700005],p_num;
18 int miu[ONE];
19 s64 Ans,sum[ONE];
20
21 int get()
22 {
23 int res=1,Q=1; char c;
24 while( (c=getchar())<48 || c>57)
25 if(c=='-')Q=-1;
26 if(Q) res=c-48;
27 while((c=getchar())>=48 && c<=57)
28 res=res*10+c-48;
29 return res*Q;
30 }
31
32 void Getmiu(int MaxN)
33 {
34 miu[1] = 1;
35 for(int i=2; i<=MaxN; i++)
36 {
37 if(!isp[i])
38 prime[++p_num] = i, miu[i] = -1;
39 for(int j=1; j<=p_num, i*prime[j]<=MaxN; j++)
40 {
41 isp[i * prime[j]] = 1;
42 if(i % prime[j] == 0)
43 {
44 miu[i * prime[j]] = 0;
45 break;
46 }
47 miu[i * prime[j]] = -miu[i];
48 }
49 }
50 for(int i=1; i<=MaxN; i++)
51 sum[i] = (sum[i-1] + (s64)i*i%MOD*(miu[i])%MOD) % MOD;
52 }
53
54 s64 Sum(int n,int m)
55 {
56 return ((s64)n*(n+1)/2%MOD) * ((s64)m*(m+1)/2%MOD) % MOD;
57 }
58
59 s64 f(int n,int m)
60 {
61 s64 Ans = 0;
62 if(n > m) swap(n,m);
63 for(int i=1, j=0; i<=n; i=j+1)
64 {
65 j = min(n/(n/i), m/(m/i));
66 Ans += Sum(n/i, m/i) * ((sum[j] - sum[i-1] + MOD) % MOD) % MOD;
67 Ans %= MOD;
68 }
69 return Ans;
70 }
71
72 int main()
73 {
74 n=get(); m=get();
75 if(n > m) swap(n,m);
76 Getmiu(n);
77 for(int i=1, j=0; i<=n; i=j+1)
78 {
79 j = min(n/(n/i), m/(m/i));
80 Ans += f(n/i, m/i) * ((s64)(i+j)*(j-i+1)/2%MOD) % MOD;
81 Ans %= MOD;
82 }
83 printf("%lld",Ans);
84 }