Luogu2260 [清华集训2012]模积和

https://www.luogu.com.cn/problem/P2260

数论分块

\[ 观察式子，(i\ne j)不好运算，我们可以通过容斥原理消除这个条件 \\ 令n\le m \\ \sum_{i=1}^n\sum_{j=1}^m (n\quad mod \quad i) \times (m\quad mod\quad i)\qquad (i\ne j) \\ =\sum_{i=1}^n (n-\lfloor \frac{n}{i} \rfloor \times i) \sum_{j=1}^m (m-\lfloor \frac{m}{j} \rfloor \times j)-\sum_{i=1}^n (n-\lfloor \frac{n}{i} \rfloor \times i)(m-\lfloor \frac{m}{i} \rfloor \times i) \\ 左边的式子简单数论分块即可，我们分析右边的式子 \\ \sum_{i=1}^n (n-\lfloor \frac{n}{i} \rfloor \times i)(m-\lfloor \frac{m}{i} \rfloor \times i) \\ =\sum_{i=n}^n (nm-m\lfloor\frac{n}{i}\rfloor i-n\lfloor\frac{m}{i}\rfloor i+\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor i^2) \\ nm,m\lfloor\frac{n}{i}\rfloor i,n\lfloor\frac{m}{i}\rfloor i 容易处理 \\ \lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor i^2 可以将n和m直接分块，与nm,m\lfloor\frac{n}{i}\rfloor i,n\lfloor\frac{m}{i}\rfloor i一起处理(具体可看代码) \\ 对于i^2，运用\sum_{i=1}^n i^2=\frac{n(n+1)(2n+1)}{6}计算 \\ 注意，19940417不是质数，可以用exgcd求出其逆元 \]

C++ Code:

#include<iostream>
#include<cstdio>
#include<cmath>
#define ll long long
#define mod 19940417
using namespace std;
int n,m,l,r;
long long inv2,inv6,ans1,ans2,ans3,ans,s1,s2,s3;
long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if (!b)
    {
        x=1;
        y=0;
        return a;
    }
    int c=exgcd(b,a%b,x,y);
    int tmp=x;
    x=y;
    y=tmp-a/b*y;
    return c;
}
long long GetN(long long w)
{
    long long x,y,z;
    z=exgcd(w,mod,x,y);
    long long t=mod/z;
    x=(x%t+t)%t;
    return x;
}
long long sum1(int x,int y)
{
    return (ll)(x+y)*(y-x+1)%mod*inv2%mod;
}
long long sum2(int x,int y)
{
    x--;
    long long sx=(ll)x*(x+1)%mod*(x*2+1)%mod*inv6%mod;
    long long sy=(ll)y*(y+1)%mod*(y*2+1)%mod*inv6%mod;
    sx=((sy-sx)%mod+mod)%mod;
    return sx;
}
long long query(int n)
{
    long long ans=(ll)n*n%mod;
    int l,r;
    for (l=1;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans-=sum1(l,r)*(n/l)%mod;
        ans=(ans%mod+mod)%mod;
    }
    return ans;
}
int main()
{
    inv2=GetN(2LL);
    inv6=GetN(6LL);
    scanf("%d%d",&n,&m);
    if (n>m)
        swap(n,m);
    ans1=query(n);
    ans2=query(m);
    ans3=(ll)n*n%mod*m%mod;
    for (l=1;l<=n;l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        s1=sum1(l,r)*(n/l)%mod*m%mod;
        s2=sum1(l,r)*(m/l)%mod*n%mod;
        s3=sum2(l,r)*(n/l)%mod*(m/l)%mod;
        ans3=ans3+s3-s1-s2;
        ans3=(ans3%mod+mod)%mod;
    }
    ans=ans1*ans2%mod;
    ans-=ans3;
    ans=(ans%mod+mod)%mod;
    cout << ans << endl;
    return 0;
}