P2260 [清华集训2012]模积和

P2260 [清华集训2012]模积和

推式子：
$$
\begin {aligned}
& \sum _ {i = 1} ^ n \sum _ {j = 1} ^ m (n \bmod i) (m \bmod j) [i \neq j] \\
= & (\sum _ {i = 1} ^ n n \bmod i) (\sum _ {j = 1} ^ m m \bmod j) - \sum _ {i = 1} ^ n (n \bmod i)(m \bmod i) \\
= & (\sum _ {i = 1} ^ n n - i \left \lfloor \frac n i \right \rfloor) (\sum _ {j = 1} ^ m m - j \left \lfloor \frac m j \right \rfloor) - \sum _ {i = 1} ^ n (n - i \left \lfloor \frac n i \right \rfloor)(m - i\left \lfloor \frac m i \right \rfloor) \\
= & (n ^ 2 - \sum _ {i = 1} ^ n i \left \lfloor \frac n i \right \rfloor) (m ^ 2 - \sum _ {j = 1} ^ m j \left \lfloor \frac m j \right \rfloor) - \sum _ {i = 1} ^ n (nm - mi \left \lfloor \frac n i \right \rfloor - ni \left \lfloor \frac m i \right \rfloor + i ^ 2 \left \lfloor \frac n i \right \rfloor \left \lfloor \frac m i \right \rfloor) \\
\end {aligned}
$$
关于计算前缀平方和，有公式：
$$
\sum _ {i = 1} ^ k = \frac {k (k + 1) (2k + 1)} 6
$$
整除分块可以在 $O(\sqrt n)$ 的时间内完成计算。

查看代码

#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long LL;
template <class Type>
void read(Type &x)
{
    char c;
    bool flag = false;
    while ((c = getchar()) < '0' || c > '9')
        c == '-' && (flag = true);
    x = c - '0';
    while ((c = getchar()) >= '0' && c <= '9')
        x = (x << 3) + (x << 1) + c - '0';
    flag && (x = ~x + 1);
}
template <class Type>
void write(Type x)
{
    x < 0 && (putchar('-'), x = ~x + 1);
    x > 9 && (write(x / 10), 0);
    putchar(x % 10 + '0');
}
const int mod = 19940417, inv2 = 9970209, inv6 = 3323403;
struct modint
{
    int k;
    modint (int _)
    {
        k = _;
    }
    friend modint operator - (modint a, modint b)
    {
        return (a.k - b.k) % mod;
    }
    friend modint operator + (modint a, modint b)
    {
        return (a.k + b.k) % mod;
    }
    friend modint operator * (modint a, modint b)
    {
        return (LL)a.k * b.k % mod;
    }
};
modint s_1 (modint k)
{
    return (k + 1) * k * inv2;
}
modint s_2 (modint k)
{
    return k * (k + 1) * (2 * k + 1) * inv6;
}
int main ()
{
    int n, m;
    read(n), read(m);
    n > m && (swap(n, m), 0);
    modint a = modint(n) * n, b = modint(m) * m, c = 0;
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        a = a - (s_1(r) - s_1(l - 1)) * (n / l);
    }
    for (int l = 1, r; l <= m; l = r + 1)
    {
        r = m / (m / l);
        b = b - (s_1(r) - s_1(l - 1)) * (m / l);
    }
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        c = c - modint(r - l + 1) * n * m + (s_1(r) - s_1(l - 1)) * (n / l) * m + (s_1(r) - s_1(l - 1)) * (m / l) * n - (s_2(r) - s_2(l - 1)) * (n / l) * (m / l);
    }
    write(((a * b + c).k + mod) % mod);
    return 0;
}