题目链接: http://www.lydsy.com/JudgeOnline/problem.php?id=2956
题目大意: 求
∑i=1n∑j=1,j≠im(nmodi)(mmodj)
∑
i
=
1
n
∑
j
=
1
,
j
≠
i
m
(
n
mod
i
)
(
m
mod
j
)
对19940417取模的值。
思路分析:
从heheda神犇的博客“heheda的数论专题练习”上翻到这样一道比较有趣的题。
一开始想歪了,一直在考虑n mod i与n mod (i+1)之间的递推关系,后来发现行不通。
既然是取模,那就可以采用“化模为除”的方法,即利用公式
nmodi=n−[n/i]∗i
n
mod
i
=
n
−
[
n
/
i
]
∗
i
化简原式,暂且不考虑
i≠j
i
≠
j
的条件可得
ans=∑i=1n∑j=1m(n−i[n/i])(m−j[m/j])=∑i=1n∑j=1m(nm−nj[m/j]−mi[n/i]+ij[n/i][m/j])=n2m2−n2∑j=1mj[m/j]−m2∑i=1ni[m/i]+∑i=1ni[n/i]∗∑j=1mj[m/j]
a
n
s
=
∑
i
=
1
n
∑
j
=
1
m
(
n
−
i
[
n
/
i
]
)
(
m
−
j
[
m
/
j
]
)
=
∑
i
=
1
n
∑
j
=
1
m
(
n
m
−
n
j
[
m
/
j
]
−
m
i
[
n
/
i
]
+
i
j
[
n
/
i
]
[
m
/
j
]
)
=
n
2
m
2
−
n
2
∑
j
=
1
m
j
[
m
/
j
]
−
m
2
∑
i
=
1
n
i
[
m
/
i
]
+
∑
i
=
1
n
i
[
n
/
i
]
∗
∑
j
=
1
m
j
[
m
/
j
]
做到这一步,既然只有i和[n/i]出现在式子里(j同理),显然可以用数论分块加速这个过程了。(至于数论分块的实现可以参见我的另一篇博客 http://blog.csdn.net/suncongbo/article/details/78819470) 分成 (2n−−√+2m−−√) ( 2 n + 2 m ) 个块,每个块内的 [n/i] [ n / i ] 与 [m/i] [ m / i ] 值均不变,块内求和直接用求和公式
getsum(n)=n(n+1)2
g
e
t
s
u
m
(
n
)
=
n
(
n
+
1
)
2
, 块内求和等于块右端点的getsum返回值减去块左端点的getsum返回值。
然后去掉 i=j i = j 的情况:
∑i=1min(n,m)(n−i[n/i])(m−i[m/i])=∑i=1min(n,m)(nm−ni[m/i]−mi[n/i]+i2[n/i][m/i])=nmmin(n,m)−∑i=1min(n,m)i(n[m/i]+m[n/i])+∑i=1min(n,m)i2[n/i][m/i]
∑
i
=
1
min
(
n
,
m
)
(
n
−
i
[
n
/
i
]
)
(
m
−
i
[
m
/
i
]
)
=
∑
i
=
1
min
(
n
,
m
)
(
n
m
−
n
i
[
m
/
i
]
−
m
i
[
n
/
i
]
+
i
2
[
n
/
i
]
[
m
/
i
]
)
=
n
m
min
(
n
,
m
)
−
∑
i
=
1
m
i
n
(
n
,
m
)
i
(
n
[
m
/
i
]
+
m
[
n
/
i
]
)
+
∑
i
=
1
m
i
n
(
n
,
m
)
i
2
[
n
/
i
]
[
m
/
i
]
, 前两项仍然可以用前面的方法进行计算。棘手的问题来了:
对于最后一项,牵扯到求连续自然数的平方和。我们知道
sqrsum(i)=∑i=1ni2=n(n+1)(2n+1)6
s
q
r
s
u
m
(
i
)
=
∑
i
=
1
n
i
2
=
n
(
n
+
1
)
(
2
n
+
1
)
6
, 但是如果把
n(n+1)(2n+1)
n
(
n
+
1
)
(
2
n
+
1
)
算出来将会很大,long long无法支持。因此只能中途对
19940417
19940417
取模。但是又要除以6,牵扯到求逆元,而19940417又不是质数,无法用费马小定理求逆元。此时做法一是通过其他方法(如exgcd)求出6对于
19940417
19940417
的逆元,二是直接算出6的逆元并作为常量使用。其实,这个数的值是
3323403
3323403
.
而我使用的是第三种做法——分类讨论。(其实是因为懒得写逆元)
观察到当 n≡0 n ≡ 0 或 n≡2(mod3) n ≡ 2 ( mod 3 ) 时, n(n+1) n ( n + 1 ) 的值可以被6整除,因此先算出 n(n+1)6 n ( n + 1 ) 6 即可。当 n≡1(mod3) n ≡ 1 ( mod 3 ) 时, 2n+1 2 n + 1 被3整除, n(n+1) n ( n + 1 ) 被2整除,故计算 2n+13 2 n + 1 3 与 n(n+1)2 n ( n + 1 ) 2 相乘即可。(当然如果模数是24而不是6我甘愿老老实实地求逆元)
注意运算符优先级问题, 不可以有任何连着三个19940417或1e9级别的数不加取模地连乘,否则必爆。
代码实现
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const int NM = 31623;
const long long P = 19940417ll;
long long bl[(NM<<2)+4];
long long tmp[(NM<<2)+4];
long long n,m;
int nm,nn,mm;
void modinc(long long &a)
{
if(a<0) a+=P;
if(a>=P) a-=P;
}
void merge(int lb1,int lb2,int rb2)
{
int i = lb1,j = lb2,k = lb1;
while(i<=lb2-1 && j<=rb2)
{
if(bl[i]<bl[j]) {tmp[k] = bl[i]; i++;}
else {tmp[k] = bl[j]; j++;}
k++;
}
while(i<=lb2-1) {tmp[k] = bl[i]; i++; k++;}
while(j<=rb2) {tmp[k] = bl[j]; j++; k++;}
for(int i=lb1; i<=rb2; i++) bl[i] = tmp[i];
}
long long sqrsum(long long rb)
{
if(rb%3==1) return ((((rb+rb+1)/3)%P)*(rb*(rb+1)/2)%P)%P;
return ((rb*(rb+1)/6)%P*(rb+rb+1))%P;
}
long long getsum(long long rb)
{
return (rb*(rb+1)/2)%P;
}
long long min_ll(long long x,long long y)
{
return x<y ? x : y;
}
int main()
{
scanf("%lld%lld",&n,&m); nm = 0;
nn = sqrt(n); mm = sqrt(m);
for(int i=1; i<=nn; i++) bl[++nm] = i;
for(int i=nn; i>=1; i--) bl[++nm] = n/i;
for(int i=1; i<=mm; i++) bl[++nm] = i;
for(int i=mm; i>=1; i--) bl[++nm] = m/i;
merge(1,1+nn+nn,nm);
long long a = 0ll,b = 0ll,c = 0ll,d = 0ll;
for(int i=1; i<=nm && bl[i]<=n; i++) {a += (getsum(bl[i])-getsum(bl[i-1])+P)*(n/bl[i])%P; modinc(a);}
for(int i=1; i<=nm && bl[i]<=m; i++) {b += (getsum(bl[i])-getsum(bl[i-1])+P)*(m/bl[i])%P; modinc(b);}
c = a*b%P; d = (min_ll(n,m)*n%P)*m%P;
for(int i=1; i<=nm && bl[i]<=min_ll(n,m); i++)
{
long long tmp = (((sqrsum(bl[i])-sqrsum(bl[i-1])+P)*(n/bl[i])%P)*(m/bl[i]))%P;
d += tmp; modinc(d);
}
for(int i=1; i<=nm && bl[i]<=min_ll(n,m); i++)
{
long long tmp = ((getsum(bl[i])-getsum(bl[i-1])+P)*((m*(n/bl[i])+n*(m/bl[i]))%P))%P;
d -= tmp; modinc(d);
}
long long ans = ((n*m)%P)*((n*m)%P)%P;
ans -= (((n*n)%P)*b)%P; modinc(ans);
ans -= (((m*m)%P)*a)%P; modinc(ans);
ans += c; modinc(ans);
ans -= d; modinc(ans);
printf("%lld\n",ans);
return 0;
}