题目传送门
我们要求
\[\sum_{i=1}^n \sum_{j=1}^m [i \ne j] \cdot (n \bmod i) \cdot (m \bmod j) \bmod 19940417
\]
注意到
\[a \bmod b=a-b\lfloor \frac ab \rfloor
\]
原式可化简为
\[\sum_{i=1}^n \sum_{j=1}^m [i \ne j] \cdot (n-i\lfloor \frac ni \rfloor) \cdot (m-j \lfloor \frac mj \rfloor)
\]
\(i \ne j\) 这个限制不好处理，我们可以容斥一下，变成
\[\sum_{i=1}^n \sum_{j=1}^m (n-i \lfloor \frac ni \rfloor) \cdot (m-j \lfloor \frac mj \rfloor)-\sum_{i=1}^{\min(n,m)} (n-i\lfloor \frac ni \rfloor) \cdot (m-i \lfloor \frac mj \rfloor)
\]
把双重 \(\sum\) 分开，变成
\[\sum_{i=1}^n (n-i \lfloor \frac ni \rfloor) \sum_{j=1}^m (m-j \lfloor \frac mj \rfloor)-\sum_{i=1}^{\min(n,m)} (n-i \lfloor \frac ni \rfloor) \cdot (m-i \lfloor \frac mi \rfloor)
\]
令
\[A=\sum_{i=1}^n (n-i \lfloor \frac ni \rfloor)
\]
\[B=\sum_{j=1}^m (m-j \lfloor \frac mj \rfloor)
\]
\[C=\sum_{i=1}^{\min(n,m)} (n-i \lfloor \frac ni \rfloor) (m-i \lfloor \frac mi \rfloor)
\]
则答案为
\[(AB-C) \bmod 19940417
\]
计算 \(A\)
把 \(n\) 提出来，变成
\[n^2-\sum_{i=1}^n i \lfloor \frac ni \rfloor
\]
后面的 \(\sum\) 整除分块即可。
计算 \(B\)
同理，把 \(m\) 提出来
\[m^2-\sum_{j=1}^m j \lfloor \frac mj \rfloor
\]
后面的 \(\sum\) 整除分块即可。
计算 \(C\)
\[\begin{align*}
C&=\sum_{i=1}^{\min(n,m)} (n-i \lfloor \frac ni \rfloor) (m-i \lfloor \frac mi \rfloor)\\
&=\sum_{i=1}^{\min(n,m)} (mn-mi \lfloor \frac ni \rfloor-ni \lfloor \frac mi \rfloor +i^2 \lfloor \frac ni \rfloor \lfloor \frac mi \rfloor)\\
&=\min(n,m) \cdot mn-m \sum_{i=1}^{\min(n,m)} i \lfloor \frac ni \rfloor-n \sum_{i=1}^{\min(n,m)} i\lfloor \frac mi \rfloor +\sum_{i=1}^{\min(n,m)} i^2 \lfloor \frac ni \rfloor \lfloor \frac mi \rfloor\\
\end{align*}
\]
二维整除分块即可。
注意到 \(19940417\) 不是质数，需要 exgcd 求逆元
code
#include <bits/stdc++.h>
#define pub public:
#define pri private:
#define fri friend:
#define Ofile(s) freopen(s".in", "r", stdin), freopen (s".out", "w", stdout)
#define Cfile(s) fclose(stdin), fclose(stdout)
#define fast ios::sync_with_stdio(false); cin.tie(NULL); cout.tie(NULL);
using namespace std;using ll = long long;
using ull = unsigned long long;
using lb = long double;
using pii = pair<int, int>;
using pll = pair<ll, ll>;
using pil = pair<int, ll>;
using pli = pair<ll, int>;constexpr int mod = 19940417;
constexpr int maxn = 2e5 + 5;ll n, m;
ll A, B, C;
ll S2, S3, S4;
ll ans;ll getmod(ll p);
ll exgcd(ll a, ll b, ll &x, ll &y);
ll inverse(ll x);
ll get_sum1(ll k);
ll get_sum2(ll k);
ll calc1(ll x, ll lim);
ll calc2(ll x, ll y, ll lim);int main() {freopen("std.in", "r", stdin);freopen("std.out", "w", stdout);fast;cin >> n >> m;A = (n % mod) * (n % mod) % mod;A = getmod(A - calc1(n, n));B = (m % mod) * (m % mod) % mod;B = getmod(B - calc1(m, m));C = (min(n, m) % mod) * (m % mod) % mod * (n % mod) % mod;S2 = calc1(m, min(n, m));S3 = calc1(n, min(n, m));S4 = calc2(n, m, min(n, m));C = getmod(C - (n % mod) * (S2 % mod) % mod);C = getmod(C - (m % mod) * (S3 % mod) % mod);C = getmod(C + S4);ans = getmod(A * B % mod - C);cout << ans;return 0;
}ll getmod(ll p){return (p % mod + mod) % mod;
}ll exgcd(ll a, ll b, ll &x, ll &y){if (!b){x = 1, y = 0;return a;} ll g = exgcd(b, a % b, y, x);y -= a / b * x;return g;
}ll inverse(ll x){ll a, b;exgcd(x, mod, a, b);return getmod(a);
}ll get_sum1(ll k){k %= mod;return k * (k + 1) % mod * inverse(2) % mod;
}ll get_sum2(ll k){k %= mod;return k * (k + 1) % mod * (2 * k + 1) % mod * inverse(6) % mod;
}ll calc1(ll x, ll lim){ll res = 0, i = 1;while (i <= lim){ll q = x / i;if (!q)break;ll j = min(x / q, lim);res = (res + (q % mod) * getmod(get_sum1(j) - get_sum1(i - 1)) % mod) % mod;i = j + 1;}return res;
}ll calc2(ll x, ll y, ll lim){ll res = 0, i = 1;while (i <= lim) {ll qx = x / i;ll qm = y / i;ll r = min(x / qx, min(y / qm, lim));res = (res + (qx % mod) * (qm % mod) % mod * getmod(get_sum2(r) - get_sum2(i - 1)) % mod) % mod;i = r + 1;}return res;
}