Description
定义 S(n,m)=∑i=1n∑j=1m(n mod i)(m mod j)[i≠j]S(n,m)=\sum\limits_{i=1}^n \sum\limits_{j=1}^m (n\bmod i)(m\bmod j)[i\ne j]S(n,m)=i=1∑nj=1∑m(nmodi)(mmodj)[i=j].
求 S(n,m) mod 19940417S(n,m)\bmod 19940417S(n,m)mod19940417 的值.
Limitations
1≤n,m≤1091\le n,m\le 10^91≤n,m≤109
1s,125MB1\text{s},125\text{MB}1s,125MB
Solution
首先,你要会两个公式:
∑i=1ni=n(n+1)2∑i=1ni2=n(n+1)(2n+1)6\sum\limits_{i=1}^n i=\frac{n(n+1)}{2}\\
\sum\limits_{i=1}^n i^2=\frac{n(n+1)(2n+1)}{6}i=1∑ni=2n(n+1)i=1∑ni2=6n(n+1)(2n+1)
下面开始推式子,假设 n≤mn\le mn≤m,否则交换 n,mn,mn,m.
先拆成两部分,提出公因式得到 (∑i=1nn mod i)×(∑j=1mm mod j)−∑i=1n(n mod i)×(m mod i)(\sum\limits_{i=1}^n n\bmod i)\times(\sum\limits_{j=1}^mm\bmod j)-\sum\limits_{i=1}^n (n\bmod i)\times (m\bmod i)(i=1∑nnmodi)×(j=1∑mmmodj)−i=1∑n(nmodi)×(mmodi)
乘号两边的式子本质相同,我们只讨论一个,注意到 ∑i=1n(n mod i)=∑i=1n(n−i⌊ni⌋)\sum\limits_{i=1}^n (n\bmod i)=\sum\limits_{i=1}^n (n-i\lfloor\frac{n}{i}\rfloor)i=1∑n(nmodi)=i=1∑n(n−i⌊in⌋)
于是可用整除分块,设当前块区间为 [l,r][l,r][l,r],则有:
∑i=lr(n mod i)=n(r−l+1)−⌊nl⌋(∑i=lri)\sum\limits_{i=l}^r (n\bmod i)=n(r-l+1)-\lfloor\frac{n}{l}\rfloor(\sum\limits_{i=l}^r i)i=l∑r(nmodi)=n(r−l+1)−⌊ln⌋(i=l∑ri)
结合开头公式,∑i=lri\sum\limits_{i=l}^r ii=l∑ri 可以 O(1)O(1)O(1) 计算,于是第一部分完成.
接下来关注第二部分,依照上面得:
∑i=1n(n mod i)×(m mod i)=∑i=1n(n−i⌊ni⌋)×(m−i⌊mi⌋)\sum_{i=1}^n (n\bmod i)\times (m\bmod i)=\sum_{i=1}^n (n-i\lfloor\frac{n}{i}\rfloor)\times (m-i\lfloor\frac{m}{i}\rfloor)i=1∑n(nmodi)×(mmodi)=i=1∑n(n−i⌊in⌋)×(m−i⌊im⌋)
把括号拆开,再重组得到:
∑i=1n(nm−i(m⌊ni⌋+n⌊mi⌋)+⌊ni⌋⌊mi⌋i2)\sum_{i=1}^n (nm-i(m\lfloor\frac{n}{i}\rfloor+n\lfloor\frac{m}{i}\rfloor)+\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor i^2)i=1∑n(nm−i(m⌊in⌋+n⌊im⌋)+⌊in⌋⌊im⌋i2)
这个式子依然能整除分块计算,同样设当前块为 [l,r][l,r][l,r]:
∑i=lr(n mod i)×(m mod i)=nm(r−l+1)+⌊nl⌋⌊ml⌋(∑i=lri2)−(m⌊nl⌋+n⌊ml⌋)(∑i=lri)\sum_{i=l}^r (n\bmod i)\times (m\bmod i)=nm(r-l+1)+\lfloor\frac{n}{l}\rfloor\lfloor\frac{m}{l}\rfloor(\sum_{i=l}^ri^2)-(m\lfloor\frac{n}{l}\rfloor+n\lfloor\frac{m}{l}\rfloor)(\sum_{i=l}^ri)i=l∑r(nmodi)×(mmodi)=nm(r−l+1)+⌊ln⌋⌊lm⌋(i=l∑ri2)−(m⌊ln⌋+n⌊lm⌋)(i=l∑ri)
不过这里每次的 rrr 计算方法不一样:
r=min(⌊n⌊nl⌋⌋,⌊m⌊ml⌋⌋)r=\min(\lfloor\frac{n}{\lfloor\frac{n}{l}\rfloor}\rfloor,\lfloor\frac{m}{\lfloor\frac{m}{l}\rfloor}\rfloor)r=min(⌊⌊ln⌋n⌋,⌊⌊lm⌋m⌋)
这样第二部分大功告成,做完了,写的时候记得步步取模.
Code
1.59KB,55ms,680KB (C++20 with O2)1.59\text{KB},55\text{ms},680\text{KB}\;\texttt{(C++20 with O2)}1.59KB,55ms,680KB(C++20 with O2)
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
using ui64 = unsigned long long;
using i128 = __int128;
using ui128 = unsigned __int128;
using f4 = float;
using f8 = double;
using f16 = long double;
template<class T>
bool chmax(T &a, const T &b){
if(a < b){ a = b; return true; }
return false;
}
template<class T>
bool chmin(T &a, const T &b){
if(a > b){ a = b; return true; }
return false;
}
constexpr int mod = 19940417, inv2 = 9970209, inv6 = 3323403;
inline int add(int x, int y) { return (x + y >= mod) ? (x + y - mod) : (x + y); }
inline int sub(int x, int y) { return add(x, mod - y); }
inline int mul(int x, int y) { return 1LL * x * y % mod; }
inline int f(int l, int r) {
int a = mul(r + 1, r), b = mul(l, l - 1);
return mul(sub(a, b), inv2);
}
inline int g(int l, int r) {
int a = mul(mul(r, r + 1), 2 * r + 1);
int b = mul(mul(l, l - 1), 2 * l - 1);
return mul(sub(a, b), inv6);
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
int n, m;
cin >> n >> m;
if (n > m) swap(n, m);
int ans = 0, prod = 0;
for (int l = 1, r; l <= m; l = r + 1) {
r = m / (m / l);
prod = add(prod, sub(mul(r - l + 1, m), mul(m / l, f(l, r))));
}
for (int l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
ans = add(ans, mul(sub(mul(r - l + 1, n), mul(n / l, f(l, r))), prod));
}
for (int l = 1, r; l <= n; l = r + 1) {
r = min(n / (n / l), m / (m / l));
int tmp = add(mul(mul(n, m), r - l + 1), mul(mul(n / l, m / l), g(l, r)));
int a = mul(m, n / l), b = mul(n, m / l);
ans = sub(ans, sub(tmp, mul(f(l, r), add(a, b))));
}
cout << ans;
return 0;
}