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

题目传送门

题意

$$\sum _{i=1}^n\sum _{j=1}^m(nmodi)(mmodj),i̸=j$$

题解

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

$$\sum _{i=1}^n\sum _{j=1}^m(nmodi)(mmodj)$$

$$=\sum _{i=1}^n(nmodi)\sum _{j=1}^m(mmodj)$$

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

$$\sum _{i=1}^n(nmodi)$$

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

$$=\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$的情况：（以下默认$n\le m$)

$$\sum _{i=1}^n(nmodi)(mmodi)$$

$$=\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 -mi\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;
}