【bzoj2956】模积和

2019-04-13 14:09发布

 求∑∑((n mod i)*(m mod j))其中1<=i<=n,1<=j<=m,i≠j。 若没有i不等于j,那么只需求sigma(n%i) * sigma(m%i)。这个可以发现商的种类只有sqrt(n)级别。枚举商发现余数是等差数列。 对于i = j,和刚才类似,考虑在n / i与m / i不变时,由两个等差数列。(s1 + i * p) * (s2 + i * q) 拆开即可。 #include #include #include #include #include #include #define Rep(i, x, y) for (int i = x; i <= y; i ++) #define RepE(i, x) for (int i = pos[x]; i; i = g[i].nex) using namespace std; typedef long long LL; const int mod = 19940417; int n, m; LL ans; LL Calc(int n) { LL ret = 0; for (int i = 1; i <= n; ) { int t = n / i, lt = n / t, st = n % lt, ed = n % i; (ret += LL(ed + st) * ((ed - st) / t + 1) / 2) %= mod; i = lt + 1; } return ret; } LL G(LL x) { LL o = x * (x + 1) / 2; if (o % 3 == 0) o = ((o / 3) % mod) * (2*x+1); else o = o % mod * (2*x + 1) / 3; return o % mod; } int main() { scanf ("%d%d", &n, &m); if (n < m) swap(n, m); ans = Calc(n) * Calc(m) % mod; for (int i = 1, lt; i <= m; i = lt + 1) { LL ta = n / i, tb = m / i; if (i <= m) lt = min(n/ta, m/tb); else lt = n / ta; LL s1 = n % lt, e1 = n % i, s2 = m % lt, e2 = m % i, l = (e1 - s1) / ta + 1; LL ret = ((l * s1 % mod) * s2 + (l*(l-1)/2 % mod) * ((ta*s2 + tb*s1) % mod) + (G(l-1) * ta % mod) * tb) % mod; ans -= ret; } ans = ((ans % mod) + mod) % mod; printf ("%lld ", ans); return 0; }