[Luogu] P2260 [清华集训2012]模积和

(Link)

Description

求(sumlimits_{i=1}^nsumlimits_{j=1}^m(nmod{i}) imes(mmod{j}),i e{j}pmod{19940417})

Solution

不妨设(nle{m})，那么

(sumlimits_{i=1}^nsumlimits_{j=1}^m(nmod{i}) imes(mmod{j}),i e{j})

(=sumlimits_{i=1}^nsumlimits_{j=1}^m(nmod{i}) imes(mmod{j})-sumlimits_{i=1}^n(nmod{i}) imes(mmod{i}))

(=sumlimits_{i=1}^n(n-ilfloorfrac{n}{i} floor)sumlimits_{j=1}^m(m-ilfloorfrac{m}{j} floor)-sumlimits_{i=1}^n(n-ilfloorfrac{n}{i} floor)(m-ilfloorfrac{m}{i} floor))

前面的一坨直接整除分块，就是后面的有点不太好搞。直接拆开：

(sumlimits_{i=1}^n(n-ilfloorfrac{n}{i} floor)(m-ilfloorfrac{m}{i} floor)

(=sumlimits_{i=1}^n(nm-ilfloorfrac{n}{i} floor{m}-ilfloorfrac{m}{i} floor{n}+i^2lfloorfrac{n}{i} floorlfloorfrac{m}{i} floor))

那不也是整除分块嘛。然后用平方和公式(sumlimits_{i=1}^ni^2=frac{n(n+1)(2n+1)}{6})就可以了。

注意模数不是质数，不能用逆元，所以要用欧拉定理求出(inv_6=3323403)。

一定要多取模，不然会爆！

Code

#include <bits/stdc++.h>

using namespace std;

#define ll long long

const ll mod = 19940417;
const ll inv6 = 3323403;

ll n, m, mul1, mul2, res;

ll read()
{
	ll x = 0ll, fl = 1ll; char ch = getchar();
	while (ch < '0' || ch > '9') { if (ch == '-') fl = -1ll; ch = getchar();}
	while (ch >= '0' && ch <= '9') {x = (x << 1ll) + (x << 3ll) + ch - '0'; ch = getchar();}
	return x * fl;
}

int main()
{
	n = read(); m = read();
	for (ll l = 1, r; l <= n; l = r + 1)
	{
		r = min(n, n / (n / l));
		mul1 = (mul1 + 1ll * n % mod * (r - l + 1ll) % mod - (n / l) % mod * ((1ll * (r - l + 1) * (l + r) / 2ll) % mod) % mod + 10ll * mod) % mod;
	}
	for (ll l = 1, r; l <= m; l = r + 1)
	{
		r = min(m, m / (m / l));
		mul2 = (mul2 + 1ll * m % mod * (r - l + 1ll) % mod - (m / l) % mod * ((1ll * (r - l + 1) * (l + r) / 2ll) % mod) % mod + 10ll * mod) % mod;
	}
	res = mul1 * mul2 % mod;
	ll d = min(n, m);
	for (ll l = 1, r; l <= d; l = r + 1)
	{
		r = min(d, min(n / (n / l), m / (m / l)));
		res = (res - 1ll * (r - l + 1) % mod * n % mod * m % mod + 1ll * (n / l) % mod * m % mod * ((1ll * (r - l + 1) * (l + r) / 2ll) % mod) % mod + 1ll * (m / l) % mod * n % mod * ((1ll * (r - l + 1) * (l + r) / 2ll) % mod) % mod - 1ll * (n / l) % mod * (m / l) % mod * (((1ll * r % mod * (r + 1) % mod * (2ll * r + 1ll) % mod * inv6 % mod - 1ll * l % mod * (l - 1) % mod * (2ll * l - 1ll) % mod * inv6 % mod + 10ll * mod) % mod) + 10ll * mod) % mod + 10ll * mod) % mod;
	}
	printf("%lld
", (res + 10ll * mod) % mod);
	return 0;
}