「BZOJ 2956」模积和
令 (l=min(n,m))。这个 (i eq j) 非常不优雅，所以我们考虑分开计算，即：
[egin{aligned}
&sum_{i=1}^{n}sum_{j=1,i
eq j}^{m}(n mod i)(mmod j)\
=&sum_{i=1}^{n}sum_{j=1}^{m}(n mod i)(mmod j)-sum_{i=1}^{ exttt{l}}(n mod i)(mmod i)\
end{aligned}
]
对于前半部分有：
[egin{aligned}
&sum_{i=1}^{n}sum_{j=1}^{m}(n mod i)(mmod j)\
=&sum_{i=1}^{n}sum_{j=1}^{m}(n -lfloor frac n i
floor*i)(m -lfloor frac m j
floor*j)\
=&sum_{i=1}^{n}(n -lfloor frac n i
floor*i)sum_{j=1}^{m}(m -lfloor frac m j
floor*j)\
=&(n^2-sum_{i=1}^{n}lfloor frac n i
floor*i)(m^2-sum_{j=1}^{m}lfloor frac m j
floor*j)\
end{aligned}
]
两次整除分块计算即可。
对于后半部分有：
[egin{aligned}
&sum_{i=1}^{l}(n mod i)(mmod i)\
=&sum_{i=1}^{l}(n -lfloor frac n i
floor*i)(m -lfloor frac m i
floor*i)\
=& lnm-nsum_{i=1}^{l}lfloor frac m i
floor*i-msum_{i=1}^{l}lfloor frac n i
floor*i+sum_{i=1}^{l}i^2lfloor frac n i
floorlfloor frac m i
floor
end{aligned}
]
根据我们小学二年级就学过的知识可以知道序列 (1^2,2^2,...,n^2) 的前缀和为 (frac{n(n+1)(2n+1)}{6})。
故对于后面三个部分每个部分做一次整除分块即可。
除以 (6) 可能需要分开除，否则可能会爆。
贴个真的很丑的代码：
/*---Author:HenryHuang---*/
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll p=19940417;
ll sum(ll n){
ll a=n,b=n+1,c=2*n+1;
if(a%2==0) a/=2;
else b/=2;
if(a%3==0) a/=3;
else if(b%3==0) b/=3;
else c/=3;
return a*b%p*c%p;
}
ll gt(ll n,ll k){
ll ans=0;
ll l,r;
for(l=1;l<=n;l=r+1){
if(l<=k) r=min(n,k/(k/l));
else r=n;
ans+=(k/l)*(l+r)*(r-l+1)/2;
}
return ans%p;
}
ll solve(ll lim,ll n,ll m){
ll ans=0;
ll l,r;
for(l=1;l<=lim;l=r+1){
if(l<=lim) r=min(lim,min(n/(n/l),m/(m/l)));
else r=lim;
(ans+=(n/l)*(m/l)%p*(sum(r)-sum(l-1)+p)%p)%p;
}
return ans%p;
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0),cout.tie(0);
ll n,m;cin>>n>>m;
ll lim=min(n,m);
ll Ans=(n*n%p-gt(n,n)+p)%p*(m*m%p-gt(m,m)+p)%p;
Ans-=lim*n%p*m%p;
Ans=(Ans+p)%p;
(Ans+=n*gt(lim,m)%p)%=p;
(Ans+=m*gt(lim,n)%p)%=p;
cout<<(Ans-solve(lim,n,m)+p)%p<<'
';
return 0;
}