[清华集训2012](洛谷2260)模积和

Posted 2020-11-03 Kananix

传送门：here

一句话题意：求$ \sum\limits_{i=1}^n\sum\limits_{j=1}^m(n\ mod\ i)(m\ mod\ j)$，$ i \neq j$

简化题意：                                                                                          

假设$ n<m$

我们令$ f(x)$表示$ \sum\limits_{i=1}^x\ x\ mod\ i$，$ f2(x,y)$表示$ \sum\limits_{i=1}^n(n\ mod\ i)(m\ mod\ i)$

则原式可化简为$ f(n)f(m)-f2(n,m)$

 

考虑求$ f(x)$                                                                                         

发现这是经典的余数求和问题 原题：here

基于数论分块,当$ i \leq x$时，$ \left\lfloor\frac{x}{i}\right\rfloor$最多只有$ 2*\sqrt{x}$种取值

因此我们枚举所有可能的值，发现在$ \left\lfloor\frac{x}{i}\right\rfloor$相同时余数是一个公差为$ \left\lfloor\frac{x}{i}\right\rfloor$的降序等差数列

直接求和即可

 

考虑求$ f2(n,m)$                                                                               

不难发现依旧是数论分块的模型

根据求$ f$的经历可以发现我们要做的是两个等差数列每一项的乘积之和

我们令数论分块的某一时刻公差$ \left\lfloor\frac{n}{i}\right\rfloor$为$ a$，$ \left\lfloor\frac{m}{i}\right\rfloor$为$ b$，项数为$ num+1$,首项分别为$ x$和$ y$

则可令第一个等差数列为$ [x]，[x-a]，[x-2*a] ... [x-num*a]$，

第二个等差数列为$ [y]，[y-b]，[y-2*b] ... [y-num*b]$

乘积结果为

$ [xy] + [xy-(ay+bx)+ab]+[xy-2(ay+bx)+4ab]+...+[xy-num(ay+bx)+num^2ab]$

即$ (num+1)*xy-\frac{num*(num+1)}{2}*(ay+bx)+\frac{num*(num+1)*(2*num+1)}{6}*ab$

需要注意的是模数不是质数不能用费马求逆元

然后就可以愉快的计算了

 

code：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define rt register int
#define ll long long
#define r read()
#define p 19940417
#define inv2 9970209
#define inv6 3323403
using namespace std;
ll read()
{
    ll x = 0; int zf = 1; char ch;
    while (ch != ‘-‘ && (ch < ‘0‘ || ch > ‘9‘)) ch = getchar();
    if (ch == ‘-‘) zf = -1, ch = getchar();
    while (ch >= ‘0‘ && ch <= ‘9‘) x = x * 10 + ch - ‘0‘, ch = getchar(); return x * zf;
}
void write(ll y)
{
    if (y < 0) putchar(‘-‘), y = -y;
    if (y > 9) write(y / 10);
    putchar(y % 10 + ‘0‘);
}
inline void writeln(ll x)
{
    write(x);putchar(‘\n‘);
}
int i,j,k,m,n,x,y,z,cnt,all,num;
ll f(const int x)//余数求和 
{
    ll ans=0;
    for(rt i=1;i<=x;)
    {
        const int R=x/(x/i);
        ans+=(ll)(x%i+x%R)*(R-i+1)>>1;
        ans%=p;i=R+1;
    }
    return ans%p;
}
ll calc(const int x,const int y,const ll a,const int b,const int num)//计算贡献 
{
    ll ans=(ll)x*y%p*num%p;
    ll sum=(ll)num*(num-1)/2;sum%=p;
    ans=ans-sum*((a*y+b*x)%p)%p;
    const ll sum2=(ll)num*(num-1)%p*(2*num-1)%p*inv6%p;
    ans=(ans+sum2*a%p*b+p)%p; 
    return ans;
}
int getv(const int x,const int y)//数论分块计算答案 
{
    ll ans=0;
    for(rt i=1;i<=x;)
    {
        const int R1=x/(x/i),R2=y/(y/i);
        const int s=min(R1,R2);
        ans+=calc(x%i,y%i,x/i,y/i,s-i+1); 
        i=s+1;
    }
    return ans%p;
}
int main()
{
    n=r;m=r;if(n>m)swap(n,m);
    writeln((f(n)*f(m)%p+p-getv(n,m))%p);
    return 0;
}