洛谷P1829 [國家集訓隊]Crash的數字表格 / JZPTAB
題意:
求$\sum_{i = 1}^{n}\sum_{j = 1}^{m}lcm(i, j)$
思路:
$\sum_{i = 1}^{n}\sum_{j = 1}^{m}lcm(i, j)$
$= \sum_{i = 1}^{n}\sum_{j = 1}^{m}\frac{i * j}{\gcd(i, j)}$
$= \sum_{d}\sum_{i = 1}^{n}\sum_{j = 1}^{m}\frac{i * j}{d}[\gcd(i, j)=d]$
$= \sum_{d}\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j = 1}^{\lfloor\frac{m}{d}\rfloor}i * j * d[\gcd(i, j)=1]$
$= \sum_{d}\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j = 1}^{\lfloor\frac{m}{d}\rfloor}i * j * d \sum_{k | gcd(i, j)} mu(k)$
$= \sum_{d} d\sum_{k} mu(k) \sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}[k|i]*i \sum_{j = 1}^{\lfloor\frac{m}{d}\rfloor}[k|j] * j$
$= \sum_{d} d\sum_{k} mu(k)\sum_{i = 1}^{\lfloor\frac{n}{dk}\rfloor}k*i \sum_{j = 1}^{\lfloor\frac{m}{dk}\rfloor}k * j$
$= \sum_{d} d\sum_{k} k^2 mu(k)\sum_{i = 1}^{\lfloor\frac{n}{dk}\rfloor} i \sum_{j = 1}^{\lfloor\frac{m}{dk}\rfloor} j$
Code:
#pragma GCC optimize(3) #pragma GCC optimize(2) #include <map> #include <set> // #include <array> #include <cmath> #include <queue> #include <stack> #include<vector> #include <cstdio> #include <cstring> #include <sstream> #include <iostream> #include <stdlib.h> #include <algorithm> // #include <unordered_map> using namespace std; typedef long long ll; typedef pair<int, int> PII; #define Time (double)clock() / CLOCKS_PER_SEC #define sd(a) scanf("%d", &a) #define sdd(a, b) scanf("%d%d", &a, &b) #define slld(a) scanf("%lld", &a) #define slldd(a, b) scanf("%lld%lld", &a, &b) const int N = 1e7 + 20; const int M = 1e5 + 20; const int mod = 20101009; const double eps = 1e-6; ll cnt = 0, primes[N], mu[N], sum[N] = {0}, sd[N] = {0}; map<ll, ll> su; bool st[N]; int n, m, p, k; void get(ll n){ mu[1] = 1; d[1] = 1; for(ll i = 2; i <= n; i ++){ if(!st[i]){ primes[cnt ++] = i; mu[i] = -1; } for(int j = 0; primes[j] <= n / i; j ++){ st[i * primes[j]] = true; if(i % primes[j] == 0){ mu[i * primes[j]] = 0; break; } mu[i * primes[j]] = - mu[i]; } } for(int i = 1; i <= n; i ++){ sum[i] = sum[i - 1] + (ll)i * i * mu[i] % mod; sum[i] %= mod; } } ll s(int n){ return (ll)n * (n + 1) / 2 % mod; } ll cal(int a, int b){ ll res = 0; if(a > b) swap(a, b); for(int l = 1, r; l <= a; l = r + 1){ r = min(a / (a / l), b / (b / l)); ll ans = 0; int ad = a / l, bd = b / l; for(int i = 1, j; i <= min(ad, bd); i = j + 1){ j = min(ad / (ad / i), bd / (bd / i)); ans += (ll)(sum[j] - sum[i - 1]) * s(ad / i) % mod * s(bd / i) % mod; ans %= mod; } res += (ll)(s(r) - s(l - 1)) * ans % mod; res %= mod; } return (res + mod) % mod; } void solve() { int n, m; sdd(n, m); printf("%lld\n", cal(n, m)); } int main() { #ifdef ONLINE_JUDGE #else freopen("/home/jungu/code/in.txt", "r", stdin); // freopen("/home/jungu/code/out.txt", "w", stdout); // freopen("/home/jungu/code/practice/out.txt","w",stdout); #endif // ios::sync_with_stdio(false); cin.tie(0), cout.tie(0); int T = 1; // sd(T); k = 10000000; get(k); // int cas = 1; while (T--) { // printf("Case #%d:", cas++); solve(); } return 0; }