数学、整除分块

求

$$\sum_{i=1}^{n} \sum_{j=1}^{m} (n \bmod i) \times (m \bmod j), i \neq j$$

$\bmod 19940417$ 的值。

感谢 gza 老师的倾情讲解。

首先 $i \neq j$ 这个条件很讨厌，去掉这个条件。

$$\sum_{i=1}^{n} \sum_{j=1}^{m} (n \bmod i) (m \bmod j) - \sum_{i=1}^{\min(n,m)} (n \bmod i) (m \bmod i)$$

钦定 $n \leq m$ 去掉 $\min$ 运算符，拆解取模运算符。

$$\sum_{i=1}^{n} \sum_{j=1}^{m} \left( n - \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \left( m - \left\lfloor \frac{m}{j} \right\rfloor \times j \right) - \sum_{i=1}^{n} \left( n - \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \left( m - \left\lfloor \frac{m}{i} \right\rfloor \times i \right)$$

处理前半部分。

$$\sum_{i=1}^{n} \sum_{j=1}^{m} \left( n - \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \left( m - \left\lfloor \frac{m}{j} \right\rfloor \times j \right)$$

$$\sum_{i=1}^{n} \sum_{j=1}^{m} nm - \sum_{i=1}^{n} \sum_{j=1}^{m} n \left( \left\lfloor \frac{m}{j} \right\rfloor \times j \right) - \sum_{i=1}^{n} \sum_{j=1}^{m} m \left( \left\lfloor \frac{n}{i} \right\rfloor \times i \right) + \sum_{i=1}^{n} \sum_{j=1}^{m} \left( \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \left( \left\lfloor \frac{m}{j} \right\rfloor \times j \right)$$

变换和化简。

$$n^2 m^2 - n^2 \sum_{j=1}^{m} \left( \left\lfloor \frac{m}{j} \right\rfloor \times j \right) - m^2 \sum_{i=1}^{n} \left( \left\lfloor \frac{n}{i} \right\rfloor \times i \right) + \sum_{i=1}^{n} \left( \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \sum_{j=1}^{m} \left( \left\lfloor \frac{m}{j} \right\rfloor \times j \right)$$

处理后半部分。

$$\sum_{i=1}^{n} \left( n - \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \left( m - \left\lfloor \frac{m}{i} \right\rfloor \times i \right)$$

$$\sum_{i=1}^{n} nm - \sum_{i=1}^{n} n \left( \left\lfloor \frac{m}{i} \right\rfloor \times i \right) - \sum_{i=1}^{n} m \left( \left\lfloor \frac{n}{i} \right\rfloor \times i \right) + \sum_{i=1}^{n} \left( \left\lfloor \frac{n}{i} \right\rfloor \times i \right) \left( \left\lfloor \frac{m}{i} \right\rfloor \times i \right)$$

$$n^2 m - n \sum_{i=1}^{n} \left( \left\lfloor \frac{m}{i} \right\rfloor \times i \right) - m \sum_{i=1}^{n} \left( \left\lfloor \frac{n}{i} \right\rfloor \times i \right) + \sum_{i=1}^{n} \left( i^2 \left\lfloor \frac{n}{i} \right\rfloor \left\lfloor \frac{m}{i} \right\rfloor \right)$$

前三部分容易计算，最后一部分需要一点手法。

由于 $\left\lfloor \dfrac{n}{i} \right\rfloor$ 和 $\left\lfloor \dfrac{m}{i} \right\rfloor$ 的取值只有 $O(\sqrt{n})$ 和 $O(\sqrt{m})$ 个，所以 $\left\lfloor \dfrac{n}{i} \right\rfloor$ 和 $\left\lfloor \dfrac{m}{i} \right\rfloor$ 的极长值相同连续段（所有 $\left\lfloor \dfrac{n}{i} \right\rfloor$ 的值相同且所有 $\left\lfloor \dfrac{m}{i} \right\rfloor$ 的值相同）数量不超过 $O(\sqrt{n} + \sqrt{m})$ 个，也可以快速计算。

由于 $n,m$ 同阶，总体时间复杂度 $O(\sqrt{n})$。

除后半式子的最后一部分外的整除分块可以用一个函数处理。

对于形如 $\sum\limits_{i=1}^{n} \left( \left\lfloor \dfrac{m}{i} \right\rfloor \times i \right)$ 的式子，传入 $n,m$ 的值即可。注意模数并非质数，平方和公式与二次方和公式中的除法需要在取模之前完成。

#include<bits/stdc++.h>
using namespace std;
constexpr int mod=19'940'417;
long long n,m,ans;
long long calculate(long long x,long long y){
    long long ret=0ll;
    for(long long l=1ll,r=1ll;l<=x&&l<=y;l=r+1){
        r=min(y/(y/l),x);
        ret+=(l+r)*(r-l+1)/2%mod*(y/l)%mod;
    }
    return ret%mod;
}
namespace pre{
    long long solution(){
        long long ret=n*n%mod*m%mod*m%mod;
        long long C2=calculate(m,m);
        long long C3=calculate(n,n);
        ret=(ret-n*n%mod*C2%mod)%mod;
        ret=(ret-m*m%mod*C3%mod)%mod;
        ret=(ret+C2*C3%mod)%mod;
        return (ret+mod)%mod;
    }
}
namespace suf{
    inline long long sum(long long x){
        static long long w[4];
        static int x2,x3;
        w[1]=x,w[2]=x+1,w[3]=2*x+1;
        for(int i=1;i<=3;i++) if(w[i]%2==0) x2=i;
        for(int i=1;i<=3;i++) if(w[i]%3==0) x3=i;
        w[x2]/=2,w[x3]/=3;
        return w[1]*w[2]%mod*w[3]%mod;
    }
    long long special_calculate(){
        long long nl=1ll,nr=1ll,ml=1ll,mr=1ll,ret=0ll;
        while(nl<=n&&ml<=n){
            nr=min(n/(n/nl),n),mr=min(m/(m/ml),n);
            ret+=(sum(min(nr,mr))-sum(max(nl,ml)-1))%mod*((n/nl)*(m/ml)%mod)%mod;
            if(nr==mr) nl=nr+1,ml=mr+1;
            else nr<mr ? nl=nr+1:ml=mr+1;
        }
        return ret%mod;
    }
    long long solution(){
        long long ret=n*n%mod*m%mod;
        ret=(ret-n*calculate(n,m)%mod)%mod;
        ret=(ret-m*calculate(n,n)%mod)%mod;
        ret=(ret+special_calculate())%mod;
        return (ret+mod)%mod;
    }
}
int main(){
    scanf("%lld%lld",&n,&m);
    if(n>m) swap(n,m);
    ans=(pre::solution()-suf::solution())%mod;
    ans=(ans+mod)%mod;
    printf("%lld\n",ans);
    return 0;
}