分类： 数论分块

题目大意：求式子：
$$\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;
}