传送门：
一句话题意：求$ \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)$
发现这是经典的余数求和问题 原题：
基于数论分块,当$ 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#include #include #include #include #define rt register int#define ll long long#define r read()#define p 19940417#define inv2 9970209#define inv6 3323403using 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;}