模积和

2019-04-13 12:19发布

题意:求求∑∑((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/im/ji * i) n/i 之类的最多只有n 个,所以分块处理
(连续自然数平方和公式请自行百度) #include #include #include using namespace std; typedef long long 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); }