「清华集训 2012」模积和「数论分块」

最新推荐文章于 2024-09-01 15:19:02 发布

_Gion

最新推荐文章于 2024-09-01 15:19:02 发布

阅读量331

点赞数

分类专栏：数学 - 数论分块

本文链接：https://blog.csdn.net/Binary_Heap/article/details/82800071

版权

数学 - 数论分块专栏收录该内容

2 篇文章 0 订阅

订阅专栏

题目传送门

题意

$\sum_{i=1}^{n}\sum_{j=1}^{m} (n \;mod \;i)(m \;mod \;j),i\ne j$

题解

先不考虑 $i = j$ 的情况。

$\sum_{i=1}^{n}\sum_{j=1}^{m} (n \;mod \;i)(m \;mod \;j)$

$=\sum_{i=1}^{n}(n \;mod \;i)\sum_{j=1}^{m} (m \;mod \;j)$

这是两个独立式子的乘积，直接考虑左边的式子，右边一样的

$\sum_{i=1}^{n}(n \;mod \;i)$

$=\sum_{i=1}^{n}(n - i\lfloor\frac{n}{i}\rfloor)$

$=n^2-\sum_{i=1}^{n}i\lfloor\frac{n}{i}\rfloor$

这就可以数论（整除）分块了.

然后考虑需要减去的 $i = j$ 的情况：（以下默认 $\leq m$ )

$\sum_{i=1}^{n} (n \;mod \;i)(m \;mod \;i)$

$=\sum_{i=1}^{n} (n - i\lfloor\frac{n}{i}\rfloor)(m - i\lfloor\frac{m}{i}\rfloor)$

$=\sum_{i=1}^{n} (nm - ni\lfloor\frac{m}{i}\rfloor - m i\lfloor\frac{n}{i}\rfloor+i^2\lfloor\frac{m}{i}\rfloor \lfloor\frac{n}{i}\rfloor)$

于是第一项直接得到，第二三项一维数论分块，第四项二维数论分块。注意 $i^2$ 的前缀和为 $\frac{n(n+1)(2n+1)}{6}$ .

#include <algorithm>
#include <cstdio>

const int MOD = 19940417;
const int inv2 = 9970209;
const int inv6 = 3323403;

int query(int n) {
	return n % MOD * 1ll * (n + 1) % MOD * inv2 % MOD;
} 

int query(int l, int r) {
	return (query(r) - query(l - 1) + MOD) % MOD;
}

int query1D(int n) {
	int ans = 0;
	for(int i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i);
		ans = (ans + 1ll * ((n / i) % MOD) % MOD * query(i, j) % MOD) % MOD;
	}
	return ans;
}

int query2(int n) {
	return n % MOD * 1ll * (n + 1) % MOD * (2ll * n + 1) % MOD * inv6 % MOD;
}

int query2(int l, int r) {
	return (query2(r) - query2(l - 1) + MOD) % MOD;
}

int query2D(int n, int m) {
	int ans = 0;
	for(int i = 1, j; i <= n; i = j + 1) {
		j = std :: min(n / (n / i), m / (m / i));
		ans = (ans + 1ll * ((n / i) % MOD) % MOD * ((m / i) % MOD) % MOD * query2(i, j) % MOD) % MOD;
	}
	return ans; 
}

int query1D_2(int n, int m) {
	int ans = 0;
	for(int i = 1, j; i <= n; i = j + 1) {
		j = m / (m / i);
		if(j > n) j = n;
		ans = (ans + 1ll * ((m / i) % MOD) % MOD * query(i, j) % MOD) % MOD;
	}
	return ans;
}

int solve(int n, int m) {
	int q1 = query1D(n), q2 = query1D(m);
	int q3 = (n % MOD * 1ll * n % MOD - q1 + MOD) % MOD;
	int q4 = (m % MOD * 1ll * m % MOD - q2 + MOD) % MOD;
	int ans = q3 * 1ll * q4 % MOD;
	int ans2 = n % MOD * 1ll * n % MOD * m % MOD;
	ans2 = (ans2 - n * 1ll * query1D_2(n, m)) % MOD;
	ans2 = (ans2 - m * 1ll * q1 % MOD) % MOD;;
	ans2 = (ans2 + query2D(n, m)) % MOD;
	ans = ((ans - ans2) % MOD + MOD) % MOD;
	return ans;
}

int main() {
	int n, m;
	scanf("%d%d", &n, &m);
	if(n > m) n ^= m ^= n ^= m; 
	printf("%d\n", solve(n, m));
	return 0;
}