【WinterCamp 2013】模积和

Description

求

$$\sum_{i=1}^{n}\sum_{j=1,j\not =i}^{m}(n\mod i)*(m\mod j)$$
n,m<=10^9，答案模19940417

Solution

首先，看到有%，心里很不爽，把它变成$n-\lfloor {n\over i}\rfloor*i$，对于$i\not =j$的情况，后面减去就可以了。我们设$n<m$
于是原式就变成了

$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n-\lfloor {n\over i}\rfloor*i)*(m-\lfloor {m\over j}\rfloor*j)-\sum_{i=1}^{n}(n-\lfloor {n\over i}\rfloor*i)*(m-\lfloor {m\over i}\rfloor*i)$$
再设 $$f(n)=\sum_{i=1}^{n}n-\lfloor {n\over i}\rfloor*i$$
于是原式就变成了 $$f(n)*f(m)-\sum_{i=1}^{n}(n-\lfloor {n\over i}\rfloor*i)*(m-\lfloor {m\over i}\rfloor*i)$$
考虑分成两部分计算

Part 1

$$f(n)=\sum_{i=1}^{n}n-\lfloor {n\over i}\rfloor*i$$
$$=n^2-\sum_{i=1}^{n}\lfloor {n\over i}\rfloor*i$$
发现 $\lfloor {n\over i}\rfloor$最多只有 $\sqrt{n}$种取值，于是我们可以使用分块。（分块大法好！）
设 $l$为一段连续的$\lfloor {n\over i}\rfloor$相等的左端点，那么这个区间的右端点 $r$便满足${n\over r}>=\lfloor {n\over l}\rfloor$，并且为最大值。移项得 $r={\lfloor {n\over \lfloor {n\over l}\rfloor}}\rfloor$。那就可以利用高斯求和搞出来了。

Part 2

$$\sum_{i=1}^{n}(n-\lfloor {n\over i}\rfloor*i)*(m-\lfloor {m\over i}\rfloor*i)$$
$$=\sum_{i=1}^{n}n*m-n*i*\lfloor {m\over i}\rfloor-m*i*\lfloor {n\over i}\rfloor+\lfloor {n\over i}\rfloor*\lfloor {m\over i}\rfloor*i^2$$
$$=n^2*m+\sum_{i=1}^{n}i*(n*\lfloor {m\over i}\rfloor+m*\lfloor {n\over i}\rfloor)-\lfloor {n\over i}\rfloor*\lfloor {m\over i}\rfloor*i^2$$
这样也就可以利用分块大法打出来了。（分块大法好！）
顺便说一下， $$\sum_{i=1}^{n}i^2=n*(n+1)*(2*n+1)\over 6$$
这里用逆元，6对于19940417的逆元是3323403。（请叫我雷锋（>_<））

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define ll long long
#define mo 19940417
#define ni 3323403
using namespace std;
ll n,m,ans;
ll sum(ll l,ll r) {return (l+r)*(r-l+1)/2%mo;}
ll mi(ll n) {
    return n*(n+1)%mo*(2*n+1)%mo*ni%mo;
}
ll getan(ll n) {
    ll l=1,r;ll an=n*n%mo;
    while (l<=n) {
        r=(n/(n/l));
        an=(an-(n/l)*sum(l,r)%mo+mo)%mo;
        l=r+1;
    }
    return an;
}
int main() {
    scanf("%lld%lld",&n,&m);
    if (n>m) swap(n,m);
    ans=(getan(n)*getan(m)%mo-n*n%mo*m%mo+mo)%mo;
    ll l=1,r;
    while (l<=n) {
        r=min(n/(n/l),m / (m / l));ll x=n / l,y=m / l;
        ans=(ans+sum(l,r)*(n*y%mo+m*x%mo)%
        mo-(mi(r)-mi(l-1)+mo)%mo*x%mo*y%mo+mo)%mo;
        l=r+1;
    }
    printf("%lld",ans);
}