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

题目传送门

题意

∑ i = 1 n ∑ j = 1 m ( n    m o d    i ) ( m    m o d    j ) , i ≠ j \sum_{i=1}^{n}\sum_{j=1}^{m} (n \;mod \;i)(m \;mod \;j),i\ne j i=1nj=1m(nmodi)(mmodj),i̸=j

题解

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

∑ i = 1 n ∑ j = 1 m ( n    m o d    i ) ( m    m o d    j ) \sum_{i=1}^{n}\sum_{j=1}^{m} (n \;mod \;i)(m \;mod \;j) i=1nj=1m(nmodi)(mmodj)

= ∑ i = 1 n ( n    m o d    i ) ∑ j = 1 m ( m    m o d    j ) =\sum_{i=1}^{n}(n \;mod \;i)\sum_{j=1}^{m} (m \;mod \;j) =i=1n(nmodi)j=1m(mmodj)

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

∑ i = 1 n ( n    m o d    i ) \sum_{i=1}^{n}(n \;mod \;i) i=1n(nmodi)

= ∑ i = 1 n ( n − i ⌊ n i ⌋ ) =\sum_{i=1}^{n}(n - i\lfloor\frac{n}{i}\rfloor) =i=1n(niin)

= ∑ i = 1 n ( n − i ⌊ n i ⌋ ) =\sum_{i=1}^{n}(n - i\lfloor\frac{n}{i}\rfloor) =i=1n(niin)

= n 2 − ∑ i = 1 n i ⌊ n i ⌋ =n^2-\sum_{i=1}^{n}i\lfloor\frac{n}{i}\rfloor =n2i=1niin

这就可以数论(整除)分块了.

然后考虑需要减去的 i = j i=j i=j的情况:(以下默认 n ≤ m n \leq m nm)

∑ i = 1 n ( n    m o d    i ) ( m    m o d    i ) \sum_{i=1}^{n} (n \;mod \;i)(m \;mod \;i) i=1n(nmodi)(mmodi)

= ∑ i = 1 n ( n − i ⌊ n i ⌋ ) ( m − i ⌊ m i ⌋ ) =\sum_{i=1}^{n} (n - i\lfloor\frac{n}{i}\rfloor)(m - i\lfloor\frac{m}{i}\rfloor) =i=1n(niin)(miim)

= ∑ i = 1 n ( n m − n i ⌊ m i ⌋ − m i ⌊ n i ⌋ + i 2 ⌊ m i ⌋ ⌊ n i ⌋ ) =\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=1n(nmniimmiin+i2imin)

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

#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;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值