题解：
$$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}}\rfloor$,$M = \lfloor{\frac{m}{d}}\rfloor$
$$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}d\sum_{t = 1}^{N}{\mu(t)} \sum_{i = 1}^{\lfloor{\frac{N}{t}}\rfloor}\sum_{j = 1}^{\lfloor{\frac{M}{t}}\rfloor}ijt^2$$
$$ans = \sum_{d = 1}^{n}d\sum_{t = 1}^{N}{\mu(t)} t^2 \sum_{i = 1}^{\lfloor{\frac{N}{t}}\rfloor}\sum_{j = 1}^{\lfloor{\frac{M}{t}}\rfloor}ij$$
$$ans = \sum_{d = 1}^{n}d\sum_{t = 1}^{N}{\mu(t)} t^2 \sum_{i = 1}^{\lfloor{\frac{N}{t}}\rfloor}i\sum_{j = 1}^{\lfloor{\frac{M}{t}}\rfloor}j$$
因为$\sum_{i = 1}^{\lfloor{\frac{N}{t}}\rfloor}i$显然是可以$O(n)$预处理的，而$\lfloor{\frac{N}{t}}\rfloor \lfloor{\frac{M}{t}}\rfloor$可以整数分块，$t^2\mu(t)$可以线性筛$O(n)$预处理，总复杂度$O(n + n\sqrt{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("\n"); 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\n", ans); 56 } 57 printf("%lld\n", ans); 58 } 59 60 int main() 61 { 62 // freopen("in.in", "r", stdin); 63 pre(); 64 work(); 65 // fclose(stdin); 66 return 0; 67 }
