∑ni=1∑mj=1[i≠j](n mod i)(m mod j)−>
∑ni=1(n−i⌊n/i⌋)∑mj=1[i≠j](m−j⌊m/j⌋)−>
∑ni=1(n−i⌊n/i⌋)∑mj=1(m−j⌊m/j⌋)−∑ni=1(n−i⌊n/i⌋)(m−i⌊m/i⌋)
其实就是把mod和i≠j化开,然后两边都用分块弄一下就好了O(n−−√)
using namespace std;
const ll Mod = 19940417;
ll n,m,N2,N6;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
if(a==0)
{
x=0,y=1;
return b;
}
ll tx,ty;
ll d=exgcd(b%a,a,tx,ty);
x=ty-b/a*tx;
y=tx;
return d;
}
ll sum1(ll x,ll y)
{
ll s1=x+y, s2=(y-x+1);
return s1%Mod*s2%Mod*N2%Mod;
}
ll sum2(ll x)
{
x%=Mod;
return x*(x+1)%Mod*(2*x+1)%Mod*N6%Mod;
}
ll sum3(ll x,ll y) { return (sum2(y)-sum2(x-1)+Mod)%Mod; }
int main()
{
ll t,d;
exgcd(2,Mod,N2,t); N2=(N2%Mod+Mod)%Mod;
d=exgcd(6,Mod,N6,t); N6=(N6%Mod+Mod)%Mod;
scanf("%lld%lld",&n,&m);
ll ret=0;
int l=1,r=m; ll s=m%Mod*m%Mod;
while(l<=r)
{
int nex=(r/(r/l));
s-=(r/l)%Mod*sum1(l,nex)%Mod;
s=(s%Mod+Mod)%Mod;
l=nex+1;
}
ret=s;
l=1,r=n; s=n%Mod*n%Mod;
while(l<=r)
{
int nex=(r/(r/l));
s-=(r/l)%Mod*sum1(l,nex)%Mod;
s=(s%Mod+Mod)%Mod;
l=nex+1;
}
ret=ret*s%Mod;
l=1,r=min(n,m); s=r%Mod*n%Mod*m%Mod;
while(l<=r)
{
int nex=min(n/(n/l),m/(m/l));
ll temp=(n/l*m+m/l*n)%Mod;
s-=sum1(l,nex)*temp%Mod;
s=(s+Mod)%Mod;
s+=sum3(l,nex)*(n/l)%Mod*(m/l)%Mod;
l=nex+1;
}
ret=((ret-s)%Mod+Mod)%Mod;
printf("%lld
",ret);
return 0;
}