[国家集训队]Crash的数字表格 / JZPTAB 莫比乌斯反演

题解：
　　$$ans = sum_{i = 1}^{n}sum_{j = 1}^{m}{frac{ij}{gcd(i, j)}}$$
　　改成枚举d(设n < m)
　　$$ans = sum_{d = 1}^{n}sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) == d]frac{ij}{d}$$
　　考虑枚举$id$
　　设$N = lfloor{frac{n}{d}} floor$,$M = lfloor{frac{m}{d}} floor$
　　$$ans = sum_{d = 1}^{n}{d}sum_{i = 1}^{N}sum_{j = 1}^{M}sum_{t | gcd(i, j)}{mu(t)ij}$$
　　把后面改成枚举$mu(t)$
　　$$ans = sum_{d = 1}^{n}dsum_{t = 1}^{N}{mu(t)} sum_{i = 1}^{lfloor{frac{N}{t}} floor}sum_{j = 1}^{lfloor{frac{M}{t}} floor}ijt^2$$
　　$$ans = sum_{d = 1}^{n}dsum_{t = 1}^{N}{mu(t)} t^2 sum_{i = 1}^{lfloor{frac{N}{t}} floor}sum_{j = 1}^{lfloor{frac{M}{t}} floor}ij$$
　　$$ans = sum_{d = 1}^{n}dsum_{t = 1}^{N}{mu(t)} t^2 sum_{i = 1}^{lfloor{frac{N}{t}} floor}isum_{j = 1}^{lfloor{frac{M}{t}} floor}j$$
　　因为$sum_{i = 1}^{lfloor{frac{N}{t}} floor}i$显然是可以$O(n)$预处理的，而$lfloor{frac{N}{t}} floor lfloor{frac{M}{t}} floor$可以整数分块，$t^2mu(t)$可以线性筛$O(n)$预处理，总复杂度$O(n + nsqrt{n})$

　　不过其实还有更优的方法，以后再说吧。。。

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 #define R register int
 4 #define AC 10000010
 5 #define p 20101009
 6 #define LL long long
 7 
 8 int n, m, N, M, tot;
 9 LL ans;
10 int prime[AC], mu[AC];
11 LL sum[AC], s[AC];//质数，mu函数，1~n求和，mu(t)*t^2求和
12 bool z[AC];
13 
14 void pre()
15 {
16     scanf("%d%d", &n, &m);
17     mu[1] = 1;
18     int b = max(n, m), x;
19     s[1] = sum[1] = 1;//初始化1，，，，因为下面的i是从2开始的
20     for(R i = 2; i <= b; i++)
21     {
22         if(!z[i]) prime[++tot] = i, mu[i] = -1;
23         s[i] = (s[i - 1] + (LL)mu[i] * ((LL)i * (LL)i)%p) % p;//error !!, mu要LL
24         sum[i] = (sum[i - 1] + i) % p;
25         for(R j = 1; j <= tot; j ++)
26         {
27             x = prime[j];
28             if(i * x > b) break;
29             z[i * x] = true;
30             if(!(i % x)) break;
31             mu[i * x] = -mu[i];
32         }
33     }
34 }
35 
36 void work()
37 {
38     if(m < n) swap(n, m);
39     //for(R i = 1; i <= m; i ++) printf("%lld ", sum[i]);
40     //printf("
");
41     LL summ;
42     for(R i = 1; i <= n; i++)//枚举d
43     {
44         int pos;
45         N = n / i, M = m / i;
46         summ = 0;
47         for(R j = 1; j <= N; j = pos + 1)//枚举t
48         {
49             pos = min(N / (N / j), M / (M / j));
50             summ += (((s[pos] - s[j - 1]) % p * sum[N / j]) % p * sum[M / j]) % p;
51             summ %= p;
52         }
53         summ = (summ * i + p) % p;//error!!!别把d给忘了！
54         ans = (ans + summ + p) % p;
55         //printf("%lld
", ans);
56     }
57     printf("%lld
", ans);
58 }
59 
60 int main()
61 {
62 //    freopen("in.in", "r", stdin);
63     pre();
64     work();
65 //    fclose(stdin);
66     return 0;
67 }

View Code