模积和

题意：求求∑∑((n mod i)*(m mod j))其中1<=i<=n,1<=j<=m,i≠j。

n,m<=10^9

分析：
不难将原式转化形式
原式=$\sum_{i=1}^{n}$(n mod i) $\sum_{j=1}^{m}$(m mod j) - $\sum_{i=1}^{min(n,m)}$(n mod i) (m mod i)
=$\sum_{i=1}^{n}$(n - i * $\lfloor n/i\rfloor$) $\sum_{j=1}^{m}$(m - j $\lfloor m/j\rfloor$) - $\sum_{i=1}^{min(n,m)}$(n - i * $\lfloor n/i\rfloor$ ) * (m - i * $\lfloor m/i\rfloor$）
=$\sum_{i=1}^{n}$(n - i * $\lfloor n/i\rfloor$) $\sum_{j=1}^{m}$(m - j $\lfloor m/j\rfloor$) - $\sum_{i=1}^{min(n,m)}$(n * m - m * i * $\lfloor n/i\rfloor$ - n * i * $\lfloor m/i\rfloor$ +$\lfloor n/i\rfloor$$\lfloor m/j\rfloor$i * i）

$\lfloor n/i\rfloor$ 之类的最多只有$\sqrt{n}$ 个，所以分块处理
(连续自然数平方和公式请自行百度)

#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;

typedef long long LL;
const LL P = 19940417,ny = 3323403;
LL n,m,ans;

LL calc(LL x,LL R) {
    LL l = 1,cnt,re = 0,r,ln,rn;
    for (;l <= R;) {
        cnt = x / l;
        r = min(x / cnt,R);
        re += cnt * (l + r) * (r - l + 1) / 2;
        l = r + 1;
    }       
    return re % P;
}

LL get(LL l,LL r) {
    LL ra = r * (r + 1) % P * (2 * r + 1) % P * ny % P;
    l --;
    LL la = l * (l + 1) % P * (2 * l + 1) % P * ny % P;
    return (ra + P - la) % P; 
}

LL work(LL R) {
    LL l = 1,nn,nm,re = 0,rn,rm,r;
    for (;l <= R;) {
        nn = n / l;
        rn = min(n / nn,R);
        nm = m / l;
        rm = min(m / nm,R);
        r = min(rn,rm);
        re = (re + nn * nm % P * get(l,r) % P) % P;
        l = r + 1;
    }       
    return re % P;
}

int main() {
    scanf("%lld%lld",&n,&m);
    LL x = (n % P * n % P + P - calc(n,n)) % P;
    LL y = (m % P * m % P + P - calc(m,m)) % P;
    ans = x * y % P;
    LL nt = min(n,m);
    ans = (ans + n % P * calc(m,nt) % P) % P;
    ans = (ans + m % P * calc(n,nt) % P) % P;
    ans = (ans + P - nt % P * n % P * m % P) % P;
    ans = (ans + P - work(nt)) % P;
    printf("%lld",ans);
}