Crash的數字表格

今天的數學課上，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只想知道表格里所有數的和mod 20101009的值。

Input

輸入的第一行包含兩個正整數，分別表示N和M。

Output

輸出一個正整數，表示表格中所有數的和mod 20101009的值。

Sample Input

4 5

Sample Output

122

【數據規模和約定】

100%的數據滿足N, M ≤ 10^7。

題解

n < m
$$\begin{array}{rl}Ans& =\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j)\\ & =\sum _{i=1}^n\sum _{j=1}^m\frac{ij}{gcd(i,j)}\\ & =\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m\frac{ij}{d}[gcd(i,j)==d]\\ & =\sum _{d=1}^n\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}ijd[gcd(i,j)==1]\\ & =\sum _{d=1}^{n}d\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}ij[gcd(i,j)==1]\\ 令:Sum(n,m)& =\sum _{i=1}^n\sum _{j=1}^{m}ij\\ f(x)& =\sum _{i=1}^n\sum _{j=1}^{m}ij[gcd(i,j)==x]\\ g(x)& =\sum _{i=1}^n\sum _{j=1}^{m}ij[x|gcd(i,j)]\\ & =x^2\sum _{i=1}^{n/x}\sum _{j=1}^{m/x}ij[1|gcd(i,j)]\\ & =\sum _{x|d}^{n}f(d)\\ 則：f(x)& =\sum _{x|d}^n\mu (\frac{d}{x})g(d)\\ f(1)& =\sum _{d=1}^n\mu (d)d^{2}Sum(\frac{n}{d},\frac{m}{d})\\ Ans& =\sum _{d=1}^{n}d\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}ij[gcd(i,j)==1]\\ & =\sum _{d=1}^{n}d\sum _{i=1}^{n/d}\mu (i)i^{2}Sum(\frac{n}{id},\frac{m}{id})\end{array}$$

• $\mu (i)i^2$用前綴和預處理
  求f(x}的時候用分塊可以在$\sqrt{n}$完成
  求Sum(x,y)的時候用分塊也可以在$\sqrt{n}$完成
  所以一次詢問可以在O(n)的時間解決
• 如果是多次詢問（10000次詢問）
  需要繼續化簡
  $$\begin{array}{rl}令T& =id，考慮對T進行分塊\\ Ans& =\sum _{d=1}^{n}d\sum _{i=1}^{n/d}\mu (i)i^{2}Sum(\frac{n}{T},\frac{m}{T})\\ & =\sum _{T=1}^{n}Sum(\frac{n}{T},\frac{m}{T})\sum _{i|T}^n\mu (i)i^2\frac{T}{i}\end{array}$$
• ${\sum }_{i|T}^n\mu (i)i^2\frac{T}{i}$是積性函數可以用線性篩預處理，令t = i × prime[j]
  • 線性篩的時候如果i % prime[j] = 0，則t對應的莫比烏斯函數為0，prime[j]的貢獻為g[i] * prime[j]
  • i % prime[j] != 0，符合積形函數直接相乘
  • i 為素數，$g[i]=i?i^2$
    Sum(x,y)可以在$\sqrt{n}$完成

// 單次詢問 O(n) #include <bits/stdc++.h> #define LL long long #define P pair<int, int> #define lowbit(x) (x & -x) #define mem(a, b) memset(a, b, sizeof(a)) #define rep(i, a, n) for (int i = a; i <= n; ++i) #define maxn 10000006 #define mid ((l + r) >> 1) #define lc rt<<1 #define rc rt<<1|1 using namespace std;int mod = 20101009; int mo[maxn], prime[maxn], g[maxn]; bool vis[maxn]; int Sum(int x) {return (1ll * (1 + x) * x / 2) % mod; } void init() {mo[1] = 1;vis[1] = 1;int len = 0;for (int i = 2; i < maxn; ++i) {if (!vis[i]) {prime[len++] = i;mo[i] = -1;}for (int j = 0; j < len && 1ll * i * prime[j] < maxn; ++j) {int t = i * prime[j];vis[t] = 1;if (i % prime[j] == 0) {mo[t] = 0;break;}mo[t] = -mo[i];}} }int solve(int n, int m) {int i = 1, j;int sum = 0;while (i <= n) {j = min(n / (n/i), m / (m/i));sum += (1ll * (g[j] - g[i-1]) * Sum(n/i) % mod) * Sum(m/i) % mod;sum %= mod;i = j + 1;}return sum; }int main() { #ifndef ONLINE_JUDGE// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout); #endifinit();for (int i = 1; i < maxn; ++i) {g[i] = 1ll * i * i * mo[i] % mod;g[i] = (g[i] + g[i-1]) % mod;}int n, m;scanf("%d %d", &n, &m);if (n > m) swap(n, m);int i = 1, j;int ans = 0;while (i <= n) {j = min(n/(n/i), m/(m/i));int t = 1ll * (i + j) * (j - i + 1) / 2 % mod;ans += 1ll * solve(n/i, m/i) * t % mod;ans %= mod;i = j + 1;}printf("%d\n", (ans + mod) % mod);return 0;} // 單次詢問O(sqrt(n)) #include <bits/stdc++.h> #define LL long long #define P pair<int, int> #define lowbit(x) (x & -x) #define mem(a, b) memset(a, b, sizeof(a)) #define rep(i, a, n) for (int i = a; i <= n; ++i) #define maxn 10000006 #define mid ((l + r) >> 1) #define lc rt<<1 #define rc rt<<1|1 using namespace std;int mod = 20101009; int mo[maxn], prime[maxn], g[maxn]; bool vis[maxn]; int Sum(int x) {return (1ll * (1 + x) * x / 2) % mod; } void init() {// mo[1] = 1;g[1] = 1;vis[1] = 1;int len = 0;for (int i = 2; i < maxn; ++i) {if (!vis[i]) {prime[len++] = i;// mo[i] = -1;g[i] = (i - 1ll * i * i % mod) % mod;}for (int j = 0; j < len && 1ll * i * prime[j] < maxn; ++j) {int t = i * prime[j];vis[i * prime[j]] = 1;if (i % prime[j] == 0) {// mo[t] = 0;g[t] = 1ll * g[i] * prime[j] % mod;break;}// mo[t] = -mo[i];g[t] = 1ll * g[i] * g[prime[j]] % mod;}}for (int i = 2; i < maxn; ++i) {g[i] = (g[i] + g[i-1]) % mod;} }int main() { #ifndef ONLINE_JUDGE// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout); #endifinit();int n, m;scanf("%d %d", &n, &m);if (n > m) swap(n, m);int i = 1, j;int ans = 0;while (i <= n) {j = min(n/(n/i), m/(m/i));ans += (1ll * Sum(n/i) * Sum(m/i) % mod) * (g[j] - g[i-1]) % mod;ans %= mod;i = j + 1;}printf("%d\n", (ans + mod) % mod);return 0; }

總結

以上是生活随笔為你收集整理的BZOJ 2154 Crash的数字表格 (莫比乌斯反演)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。