BZOJ 2154 Crash的数字表格 (莫比乌斯反演)

最新推荐文章于 2021-04-10 12:46:11 发布

原创最新推荐文章于 2021-04-10 12:46:11 发布 · 272 阅读

3 ·

CC 4.0 BY-SA版权

ACM 专栏收录该内容

260 篇文章

订阅专栏

Crash在数学课上学到最小公倍数后，尝试通过编程解决一个复杂问题：计算一个N*M的数字表格中所有数字的和，其中表格的每个格子中的数字为该行和列数的最小公倍数，结果需模20101009。

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{aligned} 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{aligned}$

$\mu(i)i^2$ 用前缀和预处理
求f(x}的时候用分块可以在 ${\sqrt{n}}$ 完成
求Sum(x,y)的时候用分块也可以在 ${\sqrt{n}}$ 完成
所以一次询问可以在O(n)的时间解决
如果是多次询问（10000次询问）
需要继续化简
$\begin{aligned} 令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{aligned}$
$\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);
#endif

    init();
    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);
#endif

    init();

    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;
}