题意:求求∑∑((n mod i)*(m mod j))其中1<=i<=n,1<=j<=m,i≠j。
n,m<=10^9
分析:
不难将原式转化形式
原式=∑ni=1(n mod i) ∑mj=1(m mod j) - ∑min(n,m)i=1(n mod i) (m mod i)
=∑ni=1(n - i * ⌊n/i⌋) ∑mj=1(m - j ⌊m/j⌋) - ∑min(n,m)i=1(n - i * ⌊n/i⌋ ) * (m - i * ⌊m/i⌋)
=∑ni=1(n - i * ⌊n/i⌋) ∑mj=1(m - j ⌊m/j⌋) - ∑min(n,m)i=1(n * m - m * i * ⌊n/i⌋ - n * i * ⌊m/i⌋ +⌊n/i⌋⌊m/j⌋i * i)
⌊n/i⌋ 之类的最多只有n−−√ 个,所以分块处理
(连续自然数平方和公式请自行百度)
#include #include #include usingnamespacestd;
typedeflonglong LL;
const LL P = 19940417,ny = 3323403;
LL n,m,ans;
LL calc(LL x,LL R) {
LL l = 1,cnt,re = 0,r,ln,rn;
for (;l <= R;) {
cnt = x / l;
r = min(x / cnt,R);
re += cnt * (l + r) * (r - l + 1) / 2;
l = r + 1;
}
return re % P;
}
LL get(LL l,LL r) {
LL ra = r * (r + 1) % P * (2 * r + 1) % P * ny % P;
l --;
LL la = l * (l + 1) % P * (2 * l + 1) % P * ny % P;
return (ra + P - la) % P;
}
LL work(LL R) {
LL l = 1,nn,nm,re = 0,rn,rm,r;
for (;l <= R;) {
nn = n / l;
rn = min(n / nn,R);
nm = m / l;
rm = min(m / nm,R);
r = min(rn,rm);
re = (re + nn * nm % P * get(l,r) % P) % P;
l = r + 1;
}
return re % P;
}
int main() {
scanf("%lld%lld",&n,&m);
LL x = (n % P * n % P + P - calc(n,n)) % P;
LL y = (m % P * m % P + P - calc(m,m)) % P;
ans = x * y % P;
LL nt = min(n,m);
ans = (ans + n % P * calc(m,nt) % P) % P;
ans = (ans + m % P * calc(n,nt) % P) % P;
ans = (ans + P - nt % P * n % P * m % P) % P;
ans = (ans + P - work(nt)) % P;
printf("%lld",ans);
}