BZOJ2956 模积和

[Solution]

The original equation equals to: sigma(n % i) * sigma(m % j) - sigma((n % i) * (m % i)). That's the solution to it. You need to divide them into sqrt(n) parts by merging each equal (n % i) and (m % i).

Another useful equation is that 1^2 + 2^2 + ... + n^2 = n * (n + 1) * (2 * n + 1) / 6.

Be careful that 19940417 is not a prime number.

[Code]

#include <cstdio>
#include <algorithm>
#include <memory.h>

using namespace std;

typedef long long qw;

#ifdef WIN32
#define lld "%I64d"
#else
#define lld "%lld"
#endif
#define _l (qw)

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

inline qw sqs(qw n) {
	return n * (n + 1) % mod * (n * 2 + 1) % mod * inv6 % mod;
}

inline qw sqn(qw n) {
	return (n + 1) * n % mod * inv2 % mod;
}

qw sovo(qw n) {
	qw s = n * n;
	for (qw i = 1, j; i <= n; i = j + 1) {
		j = n / (n / i);
		s -= (j - i + 1) * (i + j) / 2 * (n / i);
	}
	return s % mod;
}

qw sovt(qw n, qw m) {
	qw e = min(n, m);
	qw s = m * n % mod * e % mod;
	for (qw i = 1, j; i <= e; i = j + 1) {
		j = min(n / (n / i), m / (m / i));
		s -= (n * (m / i) + m * (n / i)) % mod * (sqn(j) - sqn(i - 1)) % mod;
		s += (n / i) * (m / i) % mod * (sqs(j) - sqs(i - 1) + mod) % mod;
		s %= mod; 
		if (s < 0)
			s += mod;
	}
	return s;
}

int main() {
#ifndef ONLINE_JUDGE
	freopen("in.txt", "r", stdin);
#endif

	qw n, m, ans;
	scanf(lld lld, &n, &m);
	ans = (sovo(n) * sovo(m) % mod - sovt(n, m) + mod) % mod;
	printf("%d\n", (int)ans);
}