[WC2013][JZOJ3250]模积和

题目描述

求

$$\sum_{i=1}^n\sum_{j=1\text{且}i\not=j}^m(n\ mod\ i)(m\ mod\ j)$$
其中 $1\le n,m\le 10^9$
**

题目分析

关于$i\not=j$，我们可以分开算，然后相减。
显然取模运算很烦，我们将其拆分

$$\begin{align} \sum_{i=1}^n\sum_{j=1}^m(n\ mod\ i)(m\ mod\ j)&=\sum_{i=1}^n\sum_{j=1}^m(n-\lfloor\frac{n}{i}\rfloor i)(m-\lfloor\frac{m}{j}\rfloor j)\\ &=n^2m^2-m^2\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor-n^2\sum_{i=1}^m\lfloor\frac{m}{i}\rfloor+(\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor i)(\sum_{i=1}^m\lfloor\frac{m}{i}\rfloor i) \end{align}$$
两次分块处理即可，分块怎么做参考我之前写的莫比乌斯反演学习小记的其中一个部分（ 链接在此）。
还有另一部分，令 $mi=min(n,m)$
$$\begin{align} \sum_{i=1}^{mi}(n\ mod\ i)(m\ mod\ i)&=\sum_{i=1}^{mi}(n-\lfloor\frac{n}{i}\rfloor i)(m-\lfloor\frac{m}{i}\rfloor i)\\ &=mi\times nm-m\sum_{i=1}^{mi}(\lfloor\frac{n}{i}\rfloor i)-n\sum_{i=1}^{mi}(\lfloor\frac{m}{i}\rfloor i)+\sum_{i=1}^{mi}(\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor i^2) \end{align}$$
三次分块处理。平方和公式为 $\frac{n(n+1)(2n+1)}{6}$，相信大家都能手推出来。
就是这么简单，时间复杂度 $\mathrm O(\sqrt n)$。

---

代码实现

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

typedef long long LL;

const int P=19940417;
const int six=3323403;

int n,m,mi;

int calc(int en,int x)
{
    int i=1,j,ret=0;
    while (i<=en)
    {
        j=x/(x/i);
        j=min(j,en);
        int a=((i+j)&1)?(i+j):((i+j)/2);
        int b=((j-i+1)&1)?(j-i+1):((j-i+1)/2);
        (ret+=1ll*(x/i)*a%P*b%P)%=P;
        i=j+1;
    }
    return ret;
}

int calc2()
{
    int i=1,j,ret=0;
    while (i<=mi)
    {
        int j1=n/(n/i),j2=m/(m/i);
        j=min(j1,j2);
        j=min(j,mi);
        int sum1=1ll*(j+1)*j%P*((j*2+1)%P)%P*six%P;
        int sum2=1ll*i*(i-1)%P*((i*2-1)%P)%P*six%P;
        sum1=(sum1-sum2+P)%P;
        (ret+=1ll*(n/i)*(m/i)%P*sum1%P)%=P;
        i=j+1;
    }
    return ret;
}

int main()
{
    freopen("modsum.in","r",stdin);
    freopen("modsum.out","w",stdout);
    scanf("%d %d",&n,&m);
    mi=n<m?n:m;
    int ans1=calc(n,n),ans2=calc(m,m),ans3=calc(mi,n),ans4=calc(mi,m),ans5=calc2();
    int ans=1ll*n*m%P*n%P*m%P;
    ((ans-=1ll*m*m%P*ans1%P)+=P)%=P;
    ((ans-=1ll*n*n%P*ans2%P)+=P)%=P;
    (ans+=1ll*ans1*ans2%P)%=P;
    ((ans-=1ll*mi*n%P*m%P)+=P)%=P;
    (ans+=1ll*ans3*m%P)%=P;
    (ans+=1ll*ans4*n%P)%=P;
    ((ans-=ans5)+=P)%=P;
    printf("%d\n",ans);
    fclose(stdin);
    fclose(stdout);
    return 0;
}