题意
求
$$\sum_{i=1}^n\sum_{j=1}^m (n\%i) \cdot (m\%j),i \neq j$$
$\bmod \ 19940417$的值
数据范围：$n,m \le 10^9$
题解
算法：数论分块
先不考虑条件$i \neq j$，化简原式
$$ \begin{aligned} \sum_{i=1}^n\sum_{j=1}^m (n\%i) \cdot (m\%j) &= \sum_{i=1}^n (n\%i) \cdot \sum_{j=1}^m (m\%j) \\ &= \sum_{i=1}^n (n- i\lfloor \frac{n}{i} \rfloor) \cdot \sum_{j=1}^m (m- j\lfloor \frac{m}{j} \rfloor) \\ &= (n^2-\sum_{i=1}^n i\lfloor \frac{n}{i} \rfloor) \cdot (m^2-\sum_{j=1}^m j\lfloor \frac{m}{j} \rfloor) \end{aligned} $$
对两项分别使用数论分块计算即可
inline int calc(int n){
int res=n%mod*n%mod;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
(res-=(r-l+1)*(l+r)/2%mod*(n/l)%mod)%=mod;
}
return res;
}
int main(){
ans=calc(n)*calc(m)%mod;
return 0;
}由容斥原理，现在需要减去满足$i=j$的项
$$ \begin{aligned} \sum_{i=1}^{Min(n,m)} (n\%i) \cdot (m\%i) &= \sum_{i=1}^{Min(n,m)} (n-\lfloor \frac{n}{i} \rfloor) \cdot (m-\lfloor \frac{m}{i} \rfloor) \\ &= \sum_{i=1}^{Min(n,m)} (nm-in\lfloor \frac{m}{i} \rfloor-im\lfloor \frac{n}{i} \rfloor+i^2 \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor) \\ &= nm \cdot Min(n,m)-\sum_{i=1}^{Min(n,m)}(in\lfloor \frac{m}{i} \rfloor+im\lfloor \frac{n}{i} \rfloor-i^2 \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor) \end{aligned} $$
对于求和中的前两项可以进行同样的数论分块，但对于第三项，在数论分块时会遇到累计的$i^2$，于是需要求$\sum\limits_{i=l}^r i^2$，此时需要用到公式
$$\sum_{i=1}^n i^2=\frac{n(n+1)(2n+1)}{6}$$
对于$\sum\limits_{i=l}^r i^2$，运用前缀和思想即可
而对于此公式的求值，由于取模的限制，需要先求出$6$关于模数$19940417$的逆元。此处模数为合数，不可使用费马小定理，但$(6,19940417)=1$，可以使用欧拉定理或者使用扩展欧几里得算法求得
注意：
1.数据较大，需要使用$long \ long$
2.随时取模，以免溢出
3.由于代码中有需要$/2$的地方，请将其改为乘以$2$关于模数的逆元，或将(r-l+1)*(l+r)/2写在一起以保证整除（$/6$同理，但由于分子已过大，必须先取模，因此必须使用逆元）
代码：
#include<iostream>
#include<cstdio>
using namespace std;
#define int long long
const int mod=19940417;
const int inv=3323403;
int n,m,p;
int s1,s2,s3,ans;
template <typename T>
inline T Min(T a,T b){
return a<b?a:b;
}
inline int calc(int n){
int res=n%mod*n%mod;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
(res-=(r-l+1)*(l+r)/2%mod*(n/l)%mod)%=mod;
}
return res;
}
inline int S2(int n){
return n*(n+1)%mod*(2*n+1)%mod*inv%mod;
}
signed main(){
scanf("%lld%lld",&n,&m);
ans=calc(n)*calc(m)%mod;
p=Min(n,m);
(ans-=p*n%mod*m%mod)%=mod;
for(int l=1,r;l<=p;l=r+1){
r=Min(n/(n/l),m/(m/l));
s1=(r-l+1)*(l+r)/2%mod*((n/l)*m%mod+(m/l)*n%mod)%mod;
s2=(n/l)*(m/l)%mod*(S2(r)-S2(l-1))%mod;
(ans+=s1-s2)%=mod;
}
(ans+=mod)%=mod;
printf("%lld",ans);
return 0;
}