抱歉,您的浏览器无法访问本站
本页面需要浏览器支持(启用)JavaScript
了解详情 >

今日、海を見た。もう怖くない

《退役人的自我救赎系列》——其一,也就是基础的数论知识。

整除分块

求解i=1nni\sum_{i=1}^n \frac{n}i 的值,nn101010^{10} 级别。

暴力是O(n)\mathcal{O}(n) 的,死了;但是ni\lfloor \frac{n}i\rfloor 的结果是根据块状分布的,且最多只会有2n2\sqrt{n} 种不同的值。因此,我们可以利用这个性质得到一个O(n)\mathcal{O}(\sqrt{n}) 的算法。

直接地说,ni\lfloor \frac{n}i\rfloor 的值在一段连续的区间内具有相同的值,且这个区间具有右端点nni\lfloor \frac{n}{\lfloor\frac{n}i\rfloor}\rfloor

1
2
3
4
for (int l = 1, r; l <= n; l = r + 1) {
r = min(n / (n / l), n);
procedure(l, r); // [l, r] 的值是 n / l
}

上述的写法需要注意的就是 n / l 可能为 0 导致除法除零,一般来说需要特判一下。

CQOI2007:余数求和

那来点模板题。求的是模数和,看起来和上面的整除分块没有什么关系,但是:

i=1nkmodi=i=1nkkii=nki=1nkii\sum_{i=1}^n k\mod i = \sum_{i=1}^n k - \lfloor\frac{k}i\rfloor i = nk - \sum_{i=1}^n \lfloor\frac{k}i\rfloor i

这样处理之后,右边的部分仍然不是我们熟悉的整出分块形式,所以还需要进一步处理;

我们已经知道了i=1nki\sum_{i=1}^n \lfloor\frac{k}i\rfloor 在一定范围内具有相同的值,那么在块内我们令Tj=kjT_j = \lfloor\frac{k}j\rfloor,那么就可以进行如下形式的化简:

i=1nkii=i=1nTii=ji=LRTii=LRi\sum_{i=1}^n \lfloor\frac{k}i\rfloor i = \sum_{i=1}^n T_i\cdot i = \sum_j\sum_{i=L}^R T_i\cdot \sum_{i=L}^R i

这样,在每一个分块内,左边是整除分块,右边是一个公差为 1 的等差数列;块内只需要一个定值和一个很好求的和,就可以分块来做了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
signed main() {
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
#endif
auto n = read(), k = read();
auto ans = n * k;
for (llong l = 1, r; l <= n; l = r + 1) {
if (k / l) r = min(k / (k / l), n);
else r = n;
ans -= (r - l + 1) * (l + r) / 2 * (k / l);
}
cout << ans << endl;
return 0;
}

好了,现在你已经学会了数论分块了!

清华集训2012:模积和

题目是要求了iji\ne j,所以我们需要转化为下面的形式:

i=1nj=1m[ij](n mod i)(m mod j)=i=1nj=1m(n mod i)(m mod j)i=1min(n,m)(n mod i)(m mod i)\sum_{i=1}^n\sum_{j=1}^m[i\ne j](n\ \text{mod}\ i)(m\ \text{mod}\ j) \\ = \sum_{i=1}^n\sum_{j=1}^m(n\ \text{mod}\ i)(m\ \text{mod}\ j) - \sum_{i=1}^{\min(n, m)}(n\ \text{mod}\ i)(m\ \text{mod}\ i)

然后,先处理上述式子的第一项;最直接的做法就是拆开化简:

align原式=i=1nj=1m(nnii)(mmjj)=i=1nj=1mnmniimmjjn+mjniij=n2m2nm2i=1nniin2mj=1mmjj+nmi=1nniij=1mmjj {align} \begin{aligned} 原式 &= \sum_{i=1}^n\sum_{j=1}^m (n-\lfloor\frac{n}i\rfloor i)(m - \lfloor\frac{m}j\rfloor j) \\ &= \sum_{i=1}^n\sum_{j=1}^m nm - \lfloor\frac{n}i\rfloor im - \lfloor\frac{m}j\rfloor jn + \lfloor\frac{m}j\rfloor\lfloor\frac{n}i\rfloor i j \\ &= n^2m^2 - nm^2\sum_{i=1}^n\lfloor\frac{n}i\rfloor i-n^2m\sum_{j=1}^m\lfloor\frac{m}j\rfloor j + nm\sum_{i=1}^n\lfloor\frac{n}i\rfloor i\cdot\sum_{j=1}^m\lfloor\frac{m}j\rfloor j \end{aligned}

当然,也可以类似于上面求余数和的做法那样直接拆成:

i=1n(n mod i)j=1m(m mod j)\sum_{i=1}^n(n\ \text{mod}\ i) \cdot \sum_{j=1}^m(m\ \text{mod}\ j)

也可以拆出相同的结果。

总而言之,上面的式子就被拆成了一些和上面模数求和一样的形式的组合,只需要分别求出两个互不相关的部分的和就可以算出题目种式子的左边部分;接下来处理右边部分,我们先约定k=min(n,m)k = \min(n, m)

原式=i=1k(nnii)(mmii)=i=1k(nmmniinmii+minii2)=knmmi=1kniini=1kmii+i=1knimii2\begin{aligned} 原式 &= \sum_{i=1}^k (n-\lfloor\frac{n}i\rfloor i)(m - \lfloor\frac{m}i\rfloor i) \\ &= \sum_{i=1}^k(nm - m\lfloor\frac{n}i\rfloor i - n\lfloor\frac{m}i\rfloor i + \lfloor\frac{m}i\rfloor\lfloor\frac{n}i\rfloor i^2) \\ &= knm - m\sum_{i=1}^k\lfloor\frac{n}i\rfloor i - n\sum_{i=1}^k\lfloor\frac{m}i\rfloor i + \sum_{i=1}^k\lfloor\frac{n}i\rfloor \lfloor\frac{m}i\rfloor i^2 \end{aligned}

还是一样可以使用模数分块的方法求出。

此外,提一下考研常用公式:i=1ni2=16n(n+1)(2n+1)\sum_{i=1}^n i^2 = \frac16n(n+1)(2n+1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
constexpr lll mod = 19940417;

lll fast_pow(lll a, lll b = mod - 2) {
lll ret = 1;
while (b) {
if (b & 1) (ret *= a) %= mod;
(a *= a) %= mod, b >>= 1;
}
return ret %= mod;
}

llong euler_phi(llong n) {
auto ans = n;
for (llong i = 2; i * i <= n; ++i)
if (!(n % i)) {
ans = ans / i * (i - 1);
while (!(n % i)) n /= i;
}
if (n > 1) ans = ans / n * (n - 1);
return ans;
}

signed main() {
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
#endif
lll n = read(), m = read(), k = min(n, m);
auto nm = n * m % mod, n2m = nm * n % mod;
auto n2m2 = nm * nm % mod, nm2 = nm * m % mod;
const auto calc = [](lll n, lll k) {
lll ret = 0;
minimize(n, k);
for (lll l = 1, r; l <= n; l = r + 1) {
if (k / l) r = min(k / (k / l), n);
else r = n;
auto sum = (r - l + 1) * (l + r) / 2 % mod;
(ret += sum * (k / l) % mod) %= mod;
}
return ret %= mod;
};
const auto sum_mod = [&calc](lll n, lll k) {
return (n * k % mod - calc(n, k) + mod) % mod;
};
const auto phi_mod = euler_phi(mod) - 1;
const auto inv6 = fast_pow(6, phi_mod);
const auto sum_i2 = [inv6](lll n) {
return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
};
const auto lp = sum_mod(n, n) * sum_mod(m, m) % mod;
auto knm = k * nm % mod, rp = knm;
for (lll l = 1, r; l <= k; l = r + 1) {
if (n / l && m / l)
r = min({n / (n / l), m / (m / l), k});
else r = k;
auto si2 = (sum_i2(r) - sum_i2(l - 1) + mod) % mod;
auto inm = (n / l) * (m / l) % mod;
(rp += inm * si2 % mod) %= mod;
}
auto ckn = calc(k, n), ckm = calc(k, m);
auto tmp = (m * ckn % mod + n * ckm % mod) % mod;
rp = (rp + mod - tmp) % mod;
cout << llong((lp + mod - rp) % mod) << endl;
return 0;
}

需要注意的是,19940417=7×284863119940417 = 7\times2848631,并不是一个质数,所以不能用费马小定理求出模它意义下的逆元;但是66 和它互质,故应当使用欧拉定理

欧拉定理:若gcd(a,m)=1\gcd(a, m) = 1,则aϕ(m)1 (mod m)a^{\phi(m)} \equiv 1 \ (\text{mod} \ m)

欧拉函数ϕ(n)\phi(n),表示了小于等于nn 的正整数和nn 互质的整数的个数。

标准分解式:将质因数分解的结果按照大小,由小到大排列,并将相同质因数连乘以指数形式表示。

如果一个数字的标准分解式可以写成:n=p1k1p2k2prkrn = p_1^{k_1}p_2^{k_2}\cdots p_r^{k_r},那么:

ϕ(n)=i=1rpiki1(pi1)=pnpαp1(p1)=npn(11p)\phi(n) = \prod_{i=1}^rp_i^{k_i-1}(p_i-1)=\prod_{p|n}p^{\alpha_p-1}(p-1)=n\prod_{p|n}(1-\frac1p)

可以使用上述公式计算欧拉函数。

求解单个欧拉函数值,可以使用 Pollard Rho 算法优化后根据定义求解;多个欧拉函数值则可以使用线性筛求解。

实战(?)

给了x1,x2,y1,y2x_1, x_2, y_1, y_2,求解:

i=x1x2j=y1y2(ix1+x2i+jy1+y2j)2\sum_{i=x_1}^{x_2}\sum_{j=y_1}^{y_2}(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}i\rfloor+\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}j\rfloor)^2

结果对109+710^9+7 取模。

这是 CCPC 湘潭邀请赛 2021 的 C 题,无处补题。

___MEL_9_R6YR3N0GK8937H.jpg

首先还是直接把它们乘开:

原式=(y2y1+1)i=x1x2(ix1+x2i)2+(x2x1+1)j=y1y2(jy1+y2j)2+2i=x1x2j=y1y2(ix1+x2i)(jy1+y2j)\begin{aligned} 原式 =& (y_2-y_1+1)\sum_{i=x_1}^{x_2}(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}i\rfloor)^2 + (x_2-x_1+1)\sum_{j=y_1}^{y_2}(\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}j\rfloor)^2 + \\ &2\sum_{i=x_1}^{x_2}\sum_{j=y_1}^{y_2}(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}i\rfloor)(\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}j\rfloor) \end{aligned}

可以看出前两项是一类的,第三项是另一类的;其中,前两项又可以进行这样的拆分:

原式=i=x1x2ix12+2i=x1x2ix1x2i+i=x1x2x2i2原式 = \sum_{i=x_1}^{x_2}\lfloor\frac{i}{x_1}\rfloor^2 + 2\sum_{i=x_1}^{x_2}\lfloor\frac{i}{x_1}\rfloor\lfloor\frac{x_2}i\rfloor + \sum_{i=x_1}^{x_2}\lfloor\frac{x_2}i\rfloor^2

又可以分成三项;其中第三项是我们这里提到的整除分块问题,第一项是我们再熟悉不过的“分块”,可以直接计算;那么问题就变成了中间的一部分——在这里我想要说的是这并没有什么特殊的性质(或者是我没有发现),只需要分块再分块然后直接计算就可以了。

那么回到第一组式子的第三项;注意到iijj 是互不相关的,所以直接分别计算然后乘起来就行:

原式=i=x1x2(ix1+x2i)j=y1y2(jy1+y2j)原式 = \sum_{i=x_1}^{x_2}(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}i\rfloor)\cdot\sum_{j=y_1}^{y_2}(\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}j\rfloor)

只需要分别计算后求和,之后再在模意义下乘起来就得到了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
constexpr llong mod = 1e9 + 7;

constexpr lll fast_pow(lll a, lll b = mod - 2) {
lll ret = 1;
while (b) {
if (b & 1) (ret *= a) %= mod;
(a *= a) %= mod, b >>= 1;
}
return ret %= mod;
}

constexpr lll inv6 = fast_pow(6), inv2 = fast_pow(2);

lll sum_i2(lll n) {
return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}

lll sum_i(lll n) {
return (n + 1) * n % mod * inv2 % mod;
}

lll xi2(lll l, lll n, lll k) {
lll ret = 0;
minimize(n, k);
for (lll r; l <= n; l = r + 1) {
if (k / l) r = min(k / (k / l), n);
else r = n;
auto val = (k / l) * (k / l) % mod;
(ret += (r - l + 1) * val % mod) %= mod;
}
return ret;
}

lll xi(lll l, lll n, lll k) {
lll ret = 0;
minimize(n, k);
for (lll r; l <= n; l = r + 1) {
if (k / l) r = min(k / (k / l), n);
else r = n;
(ret += (r - l + 1) * (k / l) % mod) %= mod;
}
return ret;
}

lll sum_ix(lll n, lll x) {
auto m = n % x + 1, d = n / x;
if (!d) return 0;
lll ret = x * sum_i(d - 1) % mod;
(ret += m * d % mod) %= mod;
return ret;
}

lll sum_ix2(lll n, lll x) {
auto m = n % x + 1, d = n / x;
if (!d) return 0;
lll ret = x * sum_i2(d - 1) % mod;
auto d2 = d * d % mod;
(ret += m * d2 % mod) %= mod;
return ret;
}

lll ix(lll l, lll r, lll k) {
return (sum_ix(r, k) + mod - sum_ix(l - 1, k)) % mod;
}

lll ix2(lll l, lll r, lll k) {
return (sum_ix2(r, k) + mod - sum_ix2(l - 1, k)) % mod;
}

lll ix_xi(lll l, lll n, lll x1, lll x2) {
lll ret = 0;
minimize(n, x2);
for (lll r; l <= n; l = r + 1) {
if (x2 / l) r = min(x2 / (x2 / l), n);
else r = n;
auto tmp = ix(l, r, x1);
(ret += tmp * (x2 / l) % mod) %= mod;
}
return ret;
}

signed main() {
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
#endif
for (auto T = read(); T --;) {
auto x1 = read(), x2 = read(), y1 = read(), y2 = read();
auto dy = y2 - y1 + 1, dx = x2 - x1 + 1;
auto part1 = ix2(x1, x2, x1) + xi2(x1, x2, x2) + ix_xi(x1, x2, x1, x2) * 2;
auto part2 = ix2(y1, y2, y1) + xi2(y1, y2, y2) + ix_xi(y1, y2, y1, y2) * 2;
auto part3L = (ix(x1, x2, x1) + xi(x1, x2, x2)) % mod;
auto part3R = (ix(y1, y2, y1) + xi(y1, y2, y2)) % mod;
auto part3 = part3L * part3R % mod * 2 % mod;
auto part12 = (part1 %= mod) * dy % mod + (part2 %= mod) * dx % mod;
auto ans = (part12 + part3) % mod;
cout << (llong)ans << endl;
}
return 0;
}

已通过五组测试样例。

后记

参考资料

评论