BZOJ 2956 模积和

2019-04-13 17:24发布

站内文章 / 模拟电子
15100 0 1533
题意 ： 求 ∑i=1n∑j=1m(n mod i)×(m mod j),i≠j sum_{i=1}^nsum_{j=1}^m{(nspace mod space i) imes (m space mod space j)},i eq j 最后答案对1994041719940417取模 其中 n,m≤109n,mleq10^9 sol： 最坑的地方是模数不是质数 死了。 可以把x mod yx space mod space y这玩意拆开变成x−⌊xy⌋×yx - lfloorfrac{x}{y} floor imes y 然后把整个式子拆开，先不管i≠ji eq j，可以数论分块 最后把i==ji == j的加回来。 复杂度O(n)O(sqrt{n}) #include #include #include typedef long long LL; const LL p = 19940417; LL inv2, x, y, inv6; inline LL exgcd(LL a, LL b) { if (!b) { x = 1, y = 0; return a; } LL k = exgcd(b, a % b), opt = x; x = y, y = opt - a / b * y; return k; } inline LL min(LL a, LL b) {return a > b ? b : a;} inline LL IIsum(LL l, LL r) { LL res1 = (1LL + r) % p * r % p * inv2 % p; LL res2 = (l - 1LL) % p * l % p * inv2 % p; return ((res1 - res2) % p + p) % p; } inline LL IIIsum(LL l, LL r) { LL res1 = (1LL + r) % p * r % p * (2LL * r + 1LL) % p * inv6 % p; LL res2 = l % p * (l - 1) % p * (2LL * l - 1LL) % p * inv6 % p; return ((res1 - res2) % p + p) % p; } LL n, m; int main() { exgcd(2, p), inv2 = (x % p + p) % p; exgcd(6, p), inv6 = (x % p + p) % p; scanf("%lld%lld", &n, &m); LL ans = ((n % p * n % p) * m % p * m) % p; LL opt = 0; for (LL l = 1, r; l <= n; l = r + 1) { r = n / (n / l); LL res = IIsum(l, r) * ((n / l) % p) % p; opt = (opt + res) % p; } opt = opt * m % p * m % p; ans = ((ans - opt) % p + p) % p, opt = 0; for (LL l = 1, r; l <= m; l = r + 1) { r = m / (m / l); LL res = IIsum(l, r) * ((m / l) % p) % p; opt = (opt + res) % p; } LL Msum = opt % p; opt = opt * n % p * n % p; ans = ((ans - opt) % p + p) % p, opt = 0; for (LL l = 1, r; l <= n; l = r + 1) { r = n / (n / l); LL res = IIsum(l, r) * ((n / l) % p) % p * Msum % p; opt = (opt + res) % p; }ans = ((ans + opt) % p);LL upmax = min(n, m); opt = 0; for (LL l = 1, r; l <= upmax; l = r + 1) { r = min(n / (n / l), m / (m / l)); LL res = (n / l) * (m / l) % p * IIIsum(l, r) % p; opt = (opt + res) % p; res = IIsum(l, r) * m % p * (n / l) % p; opt = ((opt - res) % p + p) % p; res = IIsum(l, r) * n % p * (m / l) % p; opt = ((opt - res) % p + p) % p; res = n * m % p * (r - l + 1LL) % p; opt = (opt + res) % p; } ans = ((ans - opt) % p + p) % p; printf ("%lld", ans); }

热门文章