题意:
求
∑i=1n∑j=1m(n mod i)∗(m mod j),i≠j
∑
i
=
1
n
∑
j
=
1
m
(
n
m
o
d
i
)
∗
(
m
m
o
d
j
)
,
i
≠
j
mod 19940417 m o d 19940417 的值
首先我们直接求是 O(n∗m) O ( n ∗ m ) 的,那么我们考虑推一下这个式子。我们发现如果没有 i≠j i ≠ j 的限制,那么我们只需要分别求出 ∑ni=1n mod i ∑ i = 1 n n m o d i 和 ∑mj=1m mod j ∑ j = 1 m m m o d j ,因为 n mod i=n−⌊ni⌋∗i n m o d i = n − ⌊ n i ⌋ ∗ i ,所以可以整数分块来算每一个的结果,可以用然后做乘积即可。但是有了这个限制我们应该怎么办呢?我们考虑用包含 i=j i = j 的总和减去 i=j i = j 的情况。
∑i=1n∑j=1m(n mod i)∗(m mod j),i≠j
∑
i
=
1
n
∑
j
=
1
m
(
n
m
o
d
i
)
∗
(
m
m
o
d
j
)
,
i
≠
j
∑i=1n∑j=1m(n mod i)∗(m mod j)−∑i=1min(n,m)(n mod i)∗(m mod i)
∑
i
=
1
n
∑
j
=
1
m
(
n
m
o
d
i
)
∗
(
m
m
o
d
j
)
−
∑
i
=
1
m
i
n
(
n
,
m
)
(
n
m
o
d
i
)
∗
(
m
m
o
d
i
)
∑i=1n(n mod i)∗∑j=1m(m mod j)−∑i=1min(n,m)(n−⌊ni⌋∗i)∗(m−⌊mi⌋∗i)
∑
i
=
1
n
(
n
m
o
d
i
)
∗
∑
j
=
1
m
(
m
m
o
d
j
)
−
∑
i
=
1
m
i
n
(
n
,
m
)
(
n
−
⌊
n
i
⌋
∗
i
)
∗
(
m
−
⌊
m
i
⌋
∗
i
)
∑i=1n(n−⌊ni⌋∗i)∗∑j=1m(m−⌊mi⌋∗i)−∑i=1min(n,m)n∗m−m∗⌊ni⌋∗i−n∗⌊mi⌋∗i+⌊ni⌋∗⌊mi⌋∗i2
∑
i
=
1
n
(
n
−
⌊
n
i
⌋
∗
i
)
∗
∑
j
=
1
m
(
m
−
⌊
m
i
⌋
∗
i
)
−
∑
i
=
1
m
i
n
(
n
,
m
)
n
∗
m
−
m
∗
⌊
n
i
⌋
∗
i
−
n
∗
⌊
m
i
⌋
∗
i
+
⌊
n
i
⌋
∗
⌊
m
i
⌋
∗
i
2
好了,推导完毕,我们发现这个式子的前半部分已经讲过可以数论分块来用 O(n−−√+m−−√) O ( n + m ) 的时间求了,后半部分同理也可以。
这里需要补充一个公式:
∑i=1ni2=16n(n+1)(2n+1)
∑
i
=
1
n
i
2
=
1
6
n
(
n
+
1
)
(
2
n
+
1
)
这里我就不给出证明了,有兴趣的话可以自行搜索。
最后吐槽一句,这种取模是真的烦,式子本来就长,再加上不停地取模操作,还不能忘了 (x+mod)%mod ( x + m o d ) % m o d ,不然还会出现负数,真的很讨厌。
代码:
#include <bits/stdc++.h>
using namespace std;
long long n,m,ans,mod=19940417,ans1,ans2,ni;
void exgcd(long long a,long long b,long long &x,long long &y)
{
if(!b)
{
x=1;
y=0;
}
else
{
exgcd(b,a%b,y,x);
y-=a/b*x;
}
}
long long calc(long long x)
{
return (x*(x+1)%mod*(2*x+1)%mod*ni%mod)%mod;
}
int main()
{
scanf("%lld%lld",&n,&m);
long long r;
long long x,y;
exgcd(6,mod,x,y);
x=(x%mod+mod)%mod;
if(!x)
x+=mod;
ni=x;
ans1=(n*n)%mod;
ans2=(m*m)%mod;
for(long long l=1;l<=n;l=r+1)
{
if(n/l!=0)
r=min(n/(n/l),n);
else
r=n;
ans1=(ans1-(((r-l+1)*(l+r)/2%mod*(n/l))%mod+mod)%mod)%mod;
ans1=(ans1+mod)%mod;
}
for(long long l=1;l<=m;l=r+1)
{
if(m/l!=0)
r=min(m/(m/l),m);
else
r=m;
ans2=(ans2-(((r-l+1)*(l+r)/2%mod*(m/l))%mod+mod)%mod)%mod;
ans2=(ans2+mod)%mod;
}
ans=(ans1*ans2)%mod;
for(int l=1;l<=min(n,m);l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans=(ans-((r-l+1)*((m*n)%mod)%mod+(calc(r)-calc(l-1))%mod*((n/l)*(m/l)%mod)%mod-((r-l+1)*(l+r)/2)%mod*m%mod*(n/l)%mod-((r-l+1)*(l+r)/2)%mod*n%mod*(m/l)%mod))%mod;
ans=(ans+mod)%mod;
}
printf("%lld\n",ans);
return 0;
}