BZOJ2956 & Luogu2260: 模积和 数论分块

Description
求∑∑((n mod i)*(m mod j))其中1<=i<=n,1<=j<=m,i≠j。

---

Sample Input
3 4

---

Sample Output
1

---

这道题数论分块，我才接触。。。
设min(n, m)为o把那个式子划一下变成
∑(n - [n / i] * i) (1<=i<=n) * ∑(m - [m / i] * i) (1<=i<=m) - ∑(n - [n / i] * i) * (m - [m / i] * i) (1<=i<=o)
(n * n - ∑([n / i] * i) (1<=i<=n)) * (m * m - ∑([m / i] * i)(1<=i<=m)) - n * m * o - ∑(n * [m / i] * i) (1<=i<=o) - ∑(m * [n / i] * i) (1<=i<=o) + ∑([n / i] * i * [m / i] * i) (1<=i<=o)
设∑([n / i] * i) (1<=i<=n)为s1，∑([m / i] * i)(1<=i<=m)为s2。
设∑(n * [m / i] * i) (1<=i<=o)为s3，∑(m * [n / i] * i) (1<=i<=o)为s4。
设∑([n / i] * [m / i] * i * i) (1<=i<=o)为s5。
i * i的前缀和这部分可以靠(i * (i + 1) * (2 * i + 1)) / 6 线性求，于是求个逆元。
不过我离线搞了，但实际上也可以不求逆元，只是我嫌麻烦。。。
s1，s2，s3，s4，s5可以数论分块搞，其他都是O(1)的。

---

#include <cstdio>
#include <cstring>

using namespace std;
typedef long long LL;
LL _min(LL x, LL y) {return x < y ? x : y;}
LL _max(LL x, LL y) {return x > y ? x : y;}
const LL mod = 19940417;
const LL phi6 = 3323403;

LL s1, s2;

int main() {
    LL n, m; scanf("%lld%lld", &n, &m);
    s1 = s2 = 0;
    for(LL l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        LL pp = (r + l) * (r - l + 1) / 2; pp %= mod;
        (s1 += (pp * (n / l)) % mod) %= mod;
    }
    for(LL l = 1, r; l <= m; l = r + 1) {
        r = m / (m / l);
        LL pp = (r + l) * (r - l + 1) / 2; pp %= mod;
        (s2 += (pp * (m / l)) % mod) %= mod;
    }
    LL nm = (n * m) % mod, uu = (nm * nm) % mod;
    (uu -= ((n * n) % mod * s2) % mod) % mod; (uu -= ((m * m) % mod * s1) % mod) % mod;
    (uu += (s1 * s2) % mod) %= mod; LL oo = _min(n, m);
    s1 = s2 = 0;
    for(LL l = 1, r; l <= oo; l = r + 1) {
        r = _min(n / (n / l), oo);
        LL pp = (r + l) * (r - l + 1) / 2; pp %= mod;
        (s1 += (pp * (n / l)) % mod) %= mod;
    }
    for(LL l = 1, r; l <= oo; l = r + 1) {
        r = _min(m / (m / l), oo);
        LL pp = (r + l) * (r - l + 1) / 2; pp %= mod;
        (s2 += (pp * (m / l)) % mod) %= mod;
    }
    (uu -= (nm * oo) % mod) % mod;
    (uu += (n * s2) % mod); (uu += (m * s1) % mod);
    LL l1 = 1, r1 = 1, l2 = 1, r2 = 1, pp = 0;
    while(1) {
        LL ll = _max(l1, l2) - 1, rr = _min(r1, r2);
        LL ls = (((ll * 2 + 1) * (ll + 1) % mod * ll) % mod * phi6) % mod;
        LL rs = (((rr * 2 + 1) * (rr + 1) % mod * rr) % mod * phi6) % mod;
        (pp += ((n / l1) * (m / l2)) % mod * (rs - ls) % mod) %= mod;
        if(r1 == oo && r2 == oo) break;
        if(r1 == r2) {
            l1 = r1 + 1;
            r1 = _min(n / (n / l1), oo);
            l2 = r2 + 1;
            r2 = _min(m / (m / l2), oo);
        }
        else if(r1 < r2) {
            l1 = r1 + 1;
            r1 = _min(n / (n / l1), oo);
        }
        else {
            l2 = r2 + 1;
            r2 = _min(m / (m / l2), oo);
        }
    }
    (uu -= pp) %= mod;
    printf("%lld\n", (uu + mod) % mod);
    return 0;
}