P2260 [清华集训2012]模积和

题意：

洛谷链接

输入 $n,m$ ，求 $\sum\limits_{i=1}^n\sum\limits_{j=1}^m(n\%i)(m\%j),i\neq j$ 。答案对 $19940417$ 取模。

看到了这道题，觉得可做，~~拿来充博客防止咕咕咕。~~

开始还在想 $19940417$ 是什么生日，结果网上搜不着，估计是出题人的生日？

然而这数字还不是个质数，我谔谔毒瘤出题人！

言归正传

由于 $(n\%i)$ 与 $(m\%j)$ 互不影响，可以单独拿出来。变为：

$Ans = \sum\limits_{i=1}^n(n\%i)\sum\limits_{j=1}^m(m\%j)$

取模不好算，转化一下就能用整除分块了： $n\%i=n-\left\lfloor\dfrac ni\right\rfloor *i$

$Ans = \sum\limits_{i=1}^n(n-\left\lfloor\frac n i\right\rfloor*i)\sum\limits_{j=1}^m(m-\left\lfloor\frac mj\right\rfloor*j)$

$Ans = (n^2-\sum\limits_{i=1}^n\left\lfloor\frac ni\right\rfloor *i)(m^2-\sum\limits_{j=1}^m\left\lfloor\frac mj\right\rfloor*j)$

测一发没过样例？然后一看漏了个条件： $i\neq j$ , 呵~

那么把 $i=j$ 的减去就可了。（默认 $n$ 小于 $m$ ）

$Jian = \sum\limits_{i=1}^n(n-\left\lfloor\frac ni\right\rfloor*i)\times(m-\left\lfloor\frac mi\right\rfloor*i)$

$Jian=\sum\limits_{i=1}^nn*m-n*\left\lfloor\frac mi\right\rfloor*i-m*\left\lfloor\frac ni\right\rfloor*i+\left\lfloor\frac ni\right\rfloor\left\lfloor\frac mi\right\rfloor*i^2$

$Jian=n*n*m-n*\sum\limits_{i=1}^n\left\lfloor\frac mi\right\rfloor *i -m*\sum\limits_{i=1}^n\left\lfloor\frac ni\right\rfloor*i+\sum\limits_{i=1}^n\left\lfloor\frac ni\right\rfloor\left\lfloor\frac mi\right\rfloor*i^2$

前三项好办，最后一项有个 $i^2$ 的前缀和，它等于 $\frac {n*(n+1)*(2n+1)} 6$ ，然后你惊奇地发现它上面爆 $long\ long$ 了，没法直接算，于是还得算出 $6$ 的逆元。

用费小结果又过不了样例了，后来一想这个毒瘤出题人模数竟然不是质数，只好打个exgcd上去。算出来6的逆元其实直接用就行了，没必要在代码里算。

感觉很长很难写？其实发现整除分块的式子都是出现好多次的，记录下答案能缩短码量。

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define QWQ cout<<"QwQ"<<endl;
#define ll long long
using namespace std;
const ll mp=19940417;

ll n,m;
ll ans1,ans2,ans3,ans4;
ll xx,yy,ni;

void exgcd(ll aa,ll bb) {
	if(!bb) { xx = 1; yy = 0; return ; }
	exgcd(bb,aa%bb);
	ll t = xx; xx = yy; yy = t - aa/bb * yy;
}

inline ll he(ll x) { return x * (x+1) %mp * (2*x+1) %mp * ni %mp; }

int main() {
	exgcd(6,mp); ni = xx;
	cin>>n>>m;
	if(n>m) swap(n,m);
	for(ll l=1,r; l<=n; l=r+1)
            r = n/(n/l), (ans1 += (n/l) * ( ((r+l)*(r-l+1)/2)%mp ) %mp) %= mp;
	for(ll l=1,r; l<=m; l=r+1)
            r = m/(m/l), (ans2 += (m/l) * ( ((r+l)*(r-l+1)/2)%mp ) %mp) %= mp;
	for(ll l=1,r; l<=n; l=r+1)
            r = min( m/(m/l), n ), (ans3 += (m/l) * ( ((r+l)*(r-l+1)/2)%mp ) %mp) %= mp;
	for(ll l=1,r; l<=n; l=r+1)
            r = min( n/(n/l), m/(m/l) ), (ans4 += (n/l) * (m/l) %mp * (he(r)-he(l-1)) %mp) %= mp;
	ll tot = (n*n-ans1) %mp * (m*m %mp-ans2) %mp;
	ll jian = ((n*n %mp *m %mp + ans4 - m*ans1 %mp - n*ans3 %mp) %mp + mp) %mp;
	cout<<(tot-jian+2*mp)%mp;
	return 0;
}