题目大意：求式子：
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j) (i \ne j)$$
的值。
那么容斥以后等价于求：
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j)-\sum_{i=1}^{min(n,m)}(n\ mod\ i)*(m\ mod\ i)$$
首先求：
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j)$$
首先，我们有定理：
$$ n\ mod\ p=n-p*\lfloor \frac{n}{p} \rfloor$$
然后配合上\(\Sigma\)的运算法则，我们可以把原式转化为：
$$n^{2}*m^{2}-m^{2}\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor-n^{2}\sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor-\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor \sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor$$
我们设：
$$ P=\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor,Q=\sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor$$
所以我们可以用数论分块来分别计算\(P,Q\)的值。
而，对于第二部分：
$$R=\sum_{i=1}^{min(n,m)}(n\ mod\ i)*(m\ mod\ i)$$
我们不妨令\(n<m\)
那么将原式展开，可得：
$$R=\sum_{i=1}^{n}n*m-n*i\lfloor \frac{m}{i} \rfloor-m*i\lfloor \frac{n}{i} \rfloor+i^{2}\lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor$$
这个式子显然也是可以数论分块的。但是我们需要计算\(i^{2}\)对于式子的贡献。
用一些简单的证明方法可以得到：
$$\sum_{i=1}^{n}i^{2}=\frac{n(2n+1)(n+1)}{6}$$
那么它的区间和即是：
$$\sum_{i=1}^{r}i^{2}-\sum_{i=1}^{l-1}i^{2}=\frac{(2r^2+2rl+2l^2-r-3l+2)(r-l+1)}{6}$$
于是我们就可以统计整个结果了。
现在我们来考虑数论分块。
数论分块是一种处理形如
$$ \sum_{i=1}^{n}\lfloor \frac{n}{i} \rfloor$$
的公式的求值的方法。
而对于函数：
$$f(x)=\lfloor \frac{a}{x} \rfloor$$
可以发现，对于连续的一段\(x\)的区间，它的值是相同的。
具体来说是，令值为\(k\)则，令区间的左端点为\(l\)，则整个区间的值为：
$$k=\lfloor \frac{n}{l} \rfloor$$
那么右端点\(r\)很显然是满足方程：
$$\lfloor \frac{n}{x} \rfloor=k$$
的最大值。
依据整除的性质得到的引理一，我们发现：
$$r=\lfloor \frac{n}{k} \rfloor$$
从而，区间为：
$$[l,\lfloor \frac{n}{\lfloor \frac{n}{l} \rfloor} \rfloor$$
附：对于引理一的证明：
引理一：求证方程：
$$\lfloor \frac{n}{x} \rfloor=k$$
的解的最大值\(t\)为
$$t=\lfloor \frac{n}{k} \rfloor$$
我们令方程的一个合法解\(x\)满足：
$$x*k+r=n\ (0\le r < x)\ (1)$$
并令：
$$t=x+d$$
那么，有：
$$\lfloor \frac{n}{x+d} \rfloor=k$$
则：
$$(x+d)*k+r’=n\ (0\le r < x+d)\ (2)$$
联立方程\((1),(2)\)可得：
$$d*k+r’=r$$
由整除的定义可得：
$$d=\lfloor \frac{r}{k} \rfloor$$
很显然，当
$$x=\lfloor \frac{n}{k} \rfloor$$
时，
$$r=n\ mod\ k,r<k,\lfloor \frac{r}{k} \rfloor=0$$
此时\(x\)取到最大值。
证毕。
另：尽量用模块化编程，而不是去化简式子。这样可以用复杂度换取正确率。
#include<iostream>
#include<cstdio>
const long long MOD = 19940417;
const long long inv2 = 9970209;
const long long inv6 = 3323403;
long long n,m;
inline long long sm1(long long A){
return (A*(A+1)%MOD)*inv2%MOD;
}
inline long long sm2(long long A){
return ((A*(A+1)%MOD)*(2*A+1)%MOD)*inv6%MOD;
}
void init(){
scanf("%lld%lld",&n,&m);
n>m?n^=m^=n^=m:0;
long long P,Q,R;
int l,r;
r=0,P=0;
while(r<n){
l=r+1;
r=(n/l)?(n/(n/l)):n;
P+=(((r-l+1)*(n/l)%MOD)*(l+r)%MOD*inv2)%MOD;
P%=MOD;
// printf("P:%lld\n",P);
}
P%=MOD;
r=0,Q=0;
while(r<m){
l=r+1;
r=(m/l)?(m/(m/l)):m;
Q+=(((r-l+1)*(m/l)%MOD)*(l+r)%MOD*inv2)%MOD;
Q%=MOD;
// printf("Q:%lld\n",Q);
}
// printf("%lld %lld\n",P,Q);
Q%=MOD;
R=(n*m%MOD)*n%MOD;
long long ans=(n*m%MOD)*(n*m%MOD)%MOD-(n*n%MOD)*Q%MOD-(m*m%MOD)*P%MOD+P*Q%MOD;
ans%=MOD;
ans+=MOD;
ans%=MOD;
r=0;
while(r<n){
l=r+1;
r=std::min((n/l)?(n/(n/l)):n,(m/l)?(m/(m/l)):m);
CNT1=m/l;
CNT2=n/l;
R-=((((sm1(r)-sm1(l-1))%MOD)*(m/l)%MOD)*n%MOD+(((sm1(r)-sm1(l-1))%MOD)*(n/l)%MOD)*m%MOD);
R+=(((sm2(r)-sm2(l-1))%MOD)*((m/l)*(n/l)%MOD))%MOD;
R%=MOD;
}
ans-=R;
ans%=MOD;
ans+=MOD;
ans%=MOD;
printf("%lld",ans);
}
int main(){
init();
return 0;
}