先算出无限制的情况，再减去i==j的情况。
无限制的情况很好算，有限制的情况需要将式子拆开。
注意最后的地方要用平方和公式，模数+1是6的倍数，于是逆元就是（模数+1）/6
#include<iostream> #include<cstring> #include<cstdlib> #include<cstdio> #define MOD(x) ((x)>=mod?(x)-mod:(x)) using namespace std; const int mod=19940417,six=3323403; int n,m,sumn,summ,l1,r1,l2,r2,l,r; void read(int &k) { int f=1;k=0;char c=getchar(); while(c<'0'||c>'9')c=='-'&&(f=-1),c=getchar(); while(c<='9'&&c>='0')k=k*10+c-'0',c=getchar(); k*=f; } int solve(int n,int m) { int sum=0; for(int i=1;i<=n;i=r+1) { l=m/(m/i+1)+1;r=m/(m/i); if(r>=n)r=n; sum=MOD(sum+(1ll*(m/i)*(r-l+1)%mod*(l+r)%mod*((mod+1)>>1)%mod)); } return sum; } int pfh(int n){return 1ll*n%mod*(n+1)%mod*(2*n+1)%mod*six%mod;} int main() { read(n);read(m); sumn=(1ll*n*n-solve(n,n))%mod;summ=(1ll*m*m-solve(m,m))%mod; int sum=1ll*min(n,m)*n%mod*m%mod; for(int i=1;i<=min(n,m);i=r+1) { r=min(n/(n/i),m/(m/i)); if(r>min(n,m))r=min(n,m); sum=MOD(sum+1ll*(n/i)*(m/i)%mod*MOD(pfh(r)+mod-pfh(i-1))%mod); } sum=(sum+mod-(1ll*m*solve(min(n,m),n)%mod)+mod-(1ll*n*solve(min(n,m),m)%mod))%mod; printf("%lld ",MOD(1ll*sumn*summ%mod+mod-sum)); }