bài viết ở blog cũ, copy lên đây ai đọc hiểu thì áp dụng
1. Giới thiệu bài toán
Nếu bạn chưa biết về bài toán Josephus thì cách tốt nhất là… Google xem nó là gì. Không khó hiểu lắm (chôm từ Wiki): Có n người đứng thành vòng tròn. Bắt đầu từ người thứ 1, đếm lần lượt tới người thứ k sẽ bị xử tử. Quy luật lập đi lập lại tới khi nào chỉ còn 1 người sống sót. Hỏi người này ở vị trí thứ mấy trong n người ban đầu?
Ví dụ ban đầu có n = 11 người, và cho k = 2, ta có 12 người đứng thành vòng tròn, bắt đầu ở người thứ 1:
--> 1
11 2
10 3
9 4
8 5
7 6
theo chiều số tăng dần, cứ đếm tới người thứ 2 thì loại: 2, 4, 6, 8, 10.
1
--> 11 X
X 3
9 X
X 5
7 X
tiếp tục, bắt đầu ở người thứ 11, loại: 1, 5, 9.
X
--> 11 .
. 3
X .
. X
7 .
tiếp tục, lần này vẫn bắt đầu ở người thứ 11, loại: 3.
.
11 .
. X
. .
. .
--> 7 .
tiếp tục, bắt đầu ở người thứ 7, loại: 11.
.
X .
. .
. .
. .
7 .
Như vậy người thứ 7 là người sống sót. Ta có W(11,2) = 7.
Daniel Erman trên Numberphile có đưa ra chi tiết cách giải cực kì hay cho k = 2: Biểu diễn n dưới dạng số nhị phân, rồi chuyển số 1 bên trái qua bên phải sẽ ra vị trí người sống sót. Ví dụ n = 11, ta có 11_{10} = 1011_2, chuyển số 1 bên trái qua bên phải ta được 0111_2 = 7_{10}. Quá lẹ, quá đẹp.
Nhưng còn với k > 2 thì sao? Có xài cách chuyển về hệ k phân rồi đảo kit đầu xuống kit cuối được không? Có lẽ là không. Trên Wiki về Josephus problem có đưa ra công thức đệ quy tổng quát, nhưng không thấy giải thích gì, thôi ta đi mò công thức vậy.
2. Hướng giải
Ta sẽ mò như trong video của Numberphile. Cho k = 3, lập bảng n, W(n) với n từ 1 tới 25 bằng code trâu:
k = 3
n | W(n)
---+-----
1 | 1
2 | 2
3 | 2
4 | 1
5 | 4
6 | 1
7 | 4
8 | 7
9 | 1
10 | 4
11 | 7
12 | 10
13 | 13
14 | 2
15 | 5
16 | 8
17 | 11
18 | 14
19 | 17
20 | 20
21 | 2
22 | 5
23 | 8
24 | 11
25 | 14
Ta thấy nó y hệt như k = 2, có điều thay vì chỉ bắt đầu với vị trí 1 rồi tăng dần theo cấp số cộng, với k = 3 chúng ta còn có thể bắt đầu ở vị trí thứ 2.
Tại sao lại lòi ra số 2 thêm ở đây? Nếu dòm lại bảng trên thì ta thấy một quy luật là
\begin{aligned}
W(n) &= W(n-1) + k\\
W(n) &= W(n) - n, &if \hspace{4pt} W(n) > n
\end{aligned}
Những giá trị n khiến W(n) > n và phải “reset” lại chinh là mấu chốt ở đây. W(n) có thể reset bằng cách trừ n, và vì chúng ta chỉ tăng thêm k vị trí mỗi lần n tăng 1, vậy với k = 3 thì W(n) sau khi reset có thể bằng 1 hoặc 2. Số 2 ở đây mà ra.
Nhưng nếu tìm theo cách đệ quy này thì độ phức tạp là O(n), trong khi với k = 2 ta có thuật toán O(\log{n}) (chuyển n thành hệ nhị phân rồi chuyển ngược về hệ thập phân). Điểm mấu chốt của thuật toán lẹ này là tìm những vị trí n mà W(n) được reset: n với W(n) < k. (lưu ý W(n) không bao giờ bằng k vì người thứ k luôn bị giết ở vòng đầu tiên)
k = 3
n | W(n)
------+--------
1 | 1
2 | 2
3 | 2
4 | 1
6 | 1
9 | 1
14 | 2
21 | 2
31 | 1
47 | 2
70 | 1
105 | 1
158 | 2
237 | 2
355 | 1
533 | 2
799 | 1
Không dễ nhìn ra được quy luật gì như k = 2. Sau đây là quá trình mò của mình:
Quy luật đầu tiên mà mình nhìn ra được, đó là n - \lfloor\frac{n}{3}\rfloor bằng với giá trị n liền trước nó, ví dụ 105 - \lfloor\frac{105}{3}\rfloor = 70, 70 - \lfloor\frac{70}{3}\rfloor = 47. Nhưng với n = 47 thì không đúng nữa. Vấn đề ở đây là khi trừ \lfloor\frac{105}{3}\rfloor là ta thực hiện vòng loại đầu tiên, rồi bắt đầu với người mới ở vị trí sau vị trí 1. Với n = 3x + 2, ta sẽ bắt đầu ở vị trí n-2. Chúng ta muốn bắt đầu ở vị trí thứ 2. Như vậy ta sẽ biểu diễn n dưới dạng 3x - 1, rồi lấy n - x, hay 47 = 3\times16 - 1, 47 - 16 = 31.
Có cái dịch ngược rồi, nhưng biết dịch ngược ở vị trí n nào? Ta phải đi từ nhỏ tới lớn, hay dịch xuôi. n có thể là 3x + 1, 3x, hay 3x - 1, trừ x đi sẽ được 2x + 1, 2x, hay 2x - 1. Như vậy với n hiện tại, ta có thể tính n liền sau bằng cách chia n cho 2 rồi nhân 3. Tới đây mình còn đi sâu vô mấy cái tiểu tiết nữa để được chính xác công thức tính n khiến W(n) reset. Nhưng nó chỉ đúng cho k = 3, với k > 3 thì lại phải mò tiếp, nên có lẽ mình sẽ không trình bày ra.
Có cách khác tổng quát hơn, mình nhận thấy qua chi tiết “chia n cho 2 rồi nhân 3”. Tức là n tăng theo cấp số nhân 1.5 lần. Tổng quát, n tăng \frac{k}{k-1} lần. Như vậy với n cho trước và n này khiến W(n) reset, ta chỉ cần xét 2 số để biết n tiếp theo, đó là n_1 = n + \lfloor\frac{n}{k-1}\rfloor và n_2 = n + \lceil\frac{n}{k-1}\rceil.
Nếu ta tìm được công thức O(1) để xét n_1 hay n_2, thì ta sẽ có công thức tổng quát với độ phức tạp xấp xỉ O(\log{n}).
- Nếu n_1 = n_2, hay (k-1) \mid n: ta khỏi cần phải xét. Và ta cũng có W(n_1) = W(n_2) = W(n).
- Nếu n_1 \neq n_2: ta đi tính W(n_1) và W(n_2), bằng công thức đệ quy ở trên, nhưng tính “ăn gian” một tí, vì ta biết các số ở giữa n và n_i cần tìm tăng theo cấp số cộng bước k:
\begin{aligned}
W(n_1) &= k(n_1 - n) + W(n)\\
W(n_2) &= k(n_2 - n) + W(n)
\end{aligned}
Rồi xét W(n_i) nào cần reset, hay W(n_i) > n_i, thì đó chính là n_i tiếp theo.
Bài toán chỉ còn lại bước tìm n khởi điểm. Bạn có thể ngạc nhiên, đương nhiên là n = 1 thì W(n) = 1 rồi, còn phải tìm khởi điểm gì nữa: vì công thức tính các “chốt trụ” phía trên chỉ đúng với n \geq k. Với n < k, công thức có tí thay đổi, vì n < k nên n_1 = n, n_2 = n + 1. Như vậy ta có các chốt từ 1 tới k-1. W(n_2) khi reset bằng cách trừ có thể phải trừ nhiều lần, vì k có thể lớn hơn n nhiều lần, nên ta sẽ xài phép chia lấy dư. Nhưng phép chia lấy dư có thể trả về 0, và ta phải reset lại thành n.
3. Code giải
C++14
using JosephusPivots = std::vector<std::pair<int,int>>;
JosephusPivots josephusPivots(int N, int k)
{
JosephusPivots p;
p.emplace_back(1, 1);
for (int i = 2; i < k && p.back().first <= N; ++i)
{
int w = (p.back().second + k) % i;
if (w == 0) w = i;
p.emplace_back(i, w);
}
while (p.back().first <= N)
{
int n = p.back().first;
int w = p.back().second;
int n1 = n + n / (k - 1);
if (n % (k - 1) == 0) { p.emplace_back(n1, w); continue; }
int n2 = n1 + 1;
int w1 = (n1 - n) * k + w;
int w2 = (n2 - n) * k + w;
if (w1 > n1) p.emplace_back(n1, w1 - n1);
else p.emplace_back(n2, w2 - n2);
}
return p;
}
Tìm được các chốt trụ rồi, ta chỉ cần binary search xem n ở đâu, nếu n có trong các chốt này thì trả về W(n) tương ứng, nếu không thì tìm chốt trụ m lớn nhất bé hơn n, rồi tính W(n) theo công thức W(n) = k(n - m) + W(m).
int josephusIndex(int n, int k, const JosephusPivots& pivots)
{
auto it = --std::upper_bound(begin(pivots), end(pivots), n,
[](int x, auto const& p){ return x < p.first; });
if (it->first == n) return it->second;
return (n - it->first) * k + it->second;
}
Full code kiểm tra tính đúng đắn với duyệt trâu:
#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <iomanip>
#include <cassert>
using JosephusPivots = std::vector<std::pair<int,int>>;
JosephusPivots josephusPivots(int, int);
int josephusIndex(int, int, const JosephusPivots&);
int main()
{
int k;
std::cout << "k = "; //for example, 18
std::cin >> k;
int N;
std::cout << "N = "; //for example, 10000
std::cin >> N;
int w = std::to_string(N).size();
std::cout << std::setw(w) << "n" << " | "
<< std::setw(w+2) << "W(n)" << "\n"
<< std::setw(w+1) << std::setfill('-') << '-' << '+'
<< std::setw(2*w+2) << "\n" << std::setfill(' ');
auto pivots = josephusPivots(N, k);
for (int n = 1; n <= N; ++n)
{
// Find W(n,k) by brute force
std::vector<int> p(n);
std::iota(begin(p), end(p), 1);
for (int i = 1; p.size() > 1; i %= k)
{
for (auto& x : p) if (i++ % k == 0) x = 0;
p.erase(std::remove(begin(p), end(p), 0), end(p));
}
// Assert W(n,k) == josephusIndex(n, k)
assert(p[0] == josephusIndex(n, k, pivots));
//
if (n % 1000 == 0)
std::cout << std::setw(w) << n << " | "
<< std::setw(w) << p[0] << " "
<< josephusIndex(n, k, pivots) << "\n";
}
}
JosephusPivots josephusPivots(int N, int k)
{
JosephusPivots p;
p.emplace_back(1, 1);
for (int i = 2; i < k && p.back().first <= N; ++i)
{
int w = (p.back().second + k) % i;
if (w == 0) w = i;
p.emplace_back(i, w);
}
while (p.back().first <= N)
{
int n = p.back().first;
int w = p.back().second;
int n1 = n + n / (k - 1);
if (n % (k - 1) == 0) { p.emplace_back(n1, w); continue; }
int n2 = n1 + 1;
int w1 = (n1 - n) * k + w;
int w2 = (n2 - n) * k + w;
if (w1 > n1) p.emplace_back(n1, w1 - n1);
else p.emplace_back(n2, w2 - n2);
}
return p;
}
int josephusIndex(int n, int k, const JosephusPivots& pivots)
{
auto it = --std::upper_bound(begin(pivots), end(pivots), n,
[](int x, auto const& p){ return x < p.first; });
if (it->first == n) return it->second;
return (n - it->first) * k + it->second;
}
4. Độ phức tạp
Độ phức tạp của thuật giải này tỉ lệ thuận với số lượng chốt trụ s mà ta tìm được trong hàm josephusPivots(n, k)
, hay là \Theta(s). Bước tìm josephusIndex
có độ phức tạp \Theta(\log{s}) nên số lượng chốt trụ s quyết định độ phức tạp.
Vì các chốt trụ này là dãy cấp số nhân: a_i = \left(\frac{k}{k-1}\right) a_{i-1} và a_1 = 1, nên ta có a_s = \left(\frac{k}{k-1}\right)^s a_1 = \left(\frac{k}{k-1}\right)^s.
Vì a_s \approx n, suy ra:
\begin{aligned}
a_s &\approx n\\
\left(\frac{k}{k-1}\right)^s &\approx n\\
s &\approx \log_{\frac{k}{k-1}}{n}\\
\end{aligned}
Vậy độ phức tạp của thuật toán là \Theta(\log_{\frac{k}{k-1}}{n})
- Với k = 2 ta có độ phức tạp là \Theta(log_2{n}) (tương đương với thuật toán trên Numberphile)
- Với k = 3 ta có độ phức tạp là \Theta(log_{1.5}{n})
- Với k = 4 ta có độ phức tạp là \Theta(log_{1.33...}{n})
- Với k = 5 ta có độ phức tạp là \Theta(log_{1.25}{n})
- v.v…
Theo công thức này thì độ phức tạp không chỉ lệ thuộc vào n mà còn lệ thuộc vào k, k càng lớn độ phức tạp càng tăng (cơ số log càng nhỏ).
Để tìm hiểu độ phức tạp khi k lớn, ta có thể viết lại công thức độ phức tạp: \Theta(\frac{\log{n}}{\log{\left(\frac{k}{k-1}\right)}})
Với k càng lớn, \log{\left(\frac{k}{k-1}\right)} càng tiến gần tới \frac{1}{k} hơn. Thật vậy, thông qua Taylor series cho \log(x) tại a = 1:
\log(x) = (x - 1) - \frac{1}{2}(x - 1)^2 + \frac{1}{3}(x - 1)^3 - \frac{1}{4}(x - 1)^4 + ...
Khi k \rightarrow +\infty, \frac{k}{k-1} \rightarrow 1+, hay nói cách khác, \left(\frac{k}{k-1} - 1\right)^i với i = 2, 3, 4, ... rất nhỏ, có thể bỏ qua. Vậy ta có:
\begin{aligned}
\log{\left(\frac{k}{k-1}\right)} &> \frac{k}{k-1} - 1\\
\log{\left(\frac{k}{k-1}\right)} &> \frac{1}{k-1} > \frac{1}{k}\\
\frac{\log{n}}{\log{\left(\frac{k}{k-1}\right)}} &< k \log{n}
\end{aligned}
Hay độ phức tạp khi k lớn là O(k\log{n}), đúng như Wiki đề cập, nhưng thiếu dẫn chứng. Ở đây ta đã chứng minh được độ phức tạp này là hoàn toàn có cơ sở.
5. Chứng minh
5.1. Công thức tính chốt trụ
UPDATE (30/12/2016)
Ta muốn chứng minh chốt trụ p_i \approx \left(\frac{k}{k-1}\right) p_{i-1}
Nếu ta chấp nhận công thức:
\begin{aligned}
W(n) =
\begin{cases}
W(n-1, k) + k & \text{if}\hspace{4pt} W(n-1, k) + k \leq n \\
W(n-1, k) + k - n & \text{otherwise}
\end{cases}
\end{aligned}
Và vị trí các chốt trụ là vị trí mà W(n) < k, gọi những giá trị này là r. Ta có:
n | W(n)
--------+--------
n | r
n + 1 | r + k
n + 2 | r + 2k
...
n + m-1| r + (m-1)k
n + m | r + mk
Ta cần tìm giá trị m sao cho:
\begin{aligned}
\begin{cases}
r + mk &> n + m\\
r + (m-1)k &\leq n + m - 1
\end{cases}
\end{aligned}
Biến đổi ta có:
\begin{aligned}
\begin{cases}
m &> \frac{n-r}{k-1}\\
m - 1 &\leq \frac{n-r}{k-1}
\end{cases}
\end{aligned}
Viết gọn lại là \frac{n-r}{k-1} < m \leq \frac{n-r}{k-1} + 1. Điều này đồng nghĩa với m = \lfloor \frac{n-r}{k-1} \rfloor + 1.
Vậy là ta ko cần “thử” n_1 hay n_2 nữa mà có thể tính thẳng luôn. Ví dụ:
- Với k = 3, n = 355, chốt trụ tiếp theo được tính thẳng: 355 + \lfloor \frac{355-1}{3-1} \rfloor + 1 = 533.
- Với k = 3, n = 533, chốt trụ tiếp theo được tính thẳng: 533 + \lfloor \frac{533-2}{3-1} \rfloor + 1 = 799.
Viết lại vòng while
của hàm josephusPivots()
:
while (p.back().first <= N)
{
int m = (p.back().first - p.back().second) / (k - 1) + 1;
p.emplace_back(p.back().first + m, p.back().second + m*k - p.back().first - m);
}
Bây giờ ta cần tìm tỉ lệ \frac{p_i}{p_{i-1}} = \frac{n+m}{n}. Ta có:
\begin{aligned}
\frac{n-r}{k-1} &< & m & &\leq & \frac{n-r}{k-1} + 1\\
n + \frac{n-r}{k-1} &< & n+m & &\leq & n + \frac{n-r}{k-1} + 1\\
1 + \frac{1-\frac{r}{n}}{k-1} &< & \frac{n+m}{n} & &\leq & 1 + \frac{1-\frac{r}{n}}{k-1} + \frac{1}{n}\\
\frac{k-\frac{r}{n}}{k-1} &< & \frac{n+m}{n} & &\leq & \frac{k-\frac{r}{n}}{k-1} + \frac{1}{n}\\
\end{aligned}
Khi n tiến tới vô cực, các phần tử có mẫu là n tiến tới 0, như vậy ta có:
\lim_{n \rightarrow \infty}{\frac{n+m}{n}} = \frac{k}{k-1}
5.2. Công thức đệ quy
UPDATE (1/1/2017)
Bây giờ ta cần chứng minh:
\begin{aligned}
W(n) =
\begin{cases}
W(n-1, k) + k & \text{if}\hspace{4pt} W(n-1, k) + k \leq n \\
W(n-1, k) + k - n & \text{otherwise}
\end{cases}
\end{aligned}
Quay lại với vòng tròn đánh số:
1
n 2
. .
. .
k+1 .
X k-1
Ta thấy sau khi người ở vị trí thứ k bị giết, bài toán tìm W(n,k) trở thành tìm W(n-1,k), nhưng thay vì bắt đầu từ 1, ta bắt đầu từ k+1. Vì vậy ta cần đổi số lại nữa:
n-k+1
n-k n-k+2
. .
2 .
1 n-2
n-1
Ta có quy luật đổi số thứ tự thứ j bên W(n-1,k) thành số thứ tự thứ i bên W(n,k) như sau:
\begin{aligned}
i =
\begin{cases}
j + k &\text{if}\hspace{4pt} j+k \leq n \\
j - n + k &\text{otherwise}
\end{cases}
\end{aligned}
Cho j = W(n-1, k), thì i = W(n,k). Thế vào ta có:
\begin{aligned}
W(n,k) =
\begin{cases}
W(n-1, k) + k &\text{if}\hspace{4pt} W(n-1, k) + k \leq n \\
W(n-1, k) + k - n &\text{otherwise}
\end{cases}
\end{aligned}
Đúng y hệt công thức trên. Như vậy cách giải của chúng ta là đúng đắn.
5.3. Công thức mò
UPDATE (5/1/2017)
Chứng minh cách mò n_1 = n + \lfloor\frac{n}{k-1}\rfloor và n_2 = n + \lceil\frac{n}{k-1}\rceil.
Ta có:
\begin{aligned}
\frac{n-r}{k-1} &< & m &&\leq & \frac{n-r}{k-1} + 1\\
\frac{n}{k-1} - \frac{r}{k-1} &< & m &&\leq & \frac{n}{k-1} - \frac{r}{k-1} + 1\\
\end{aligned}
Ta đặt \frac{n}{k-1} = \lfloor\frac{n}{k-1}\rfloor + 0.R với 0.R là phần số thập phân, 0 \leq 0.R < 1.
Vì 0 < r < k nên 0 < \frac{r}{k-1} \leq 1. Đặt 0.X = 0.R - \frac{r}{k-1}, ta có -1 \leq 0.X < 1.
Vậy ta có \lfloor\frac{n}{k-1}\rfloor + 0.X < m \leq \lfloor\frac{n}{k-1}\rfloor + 0.X + 1 hay viết tắt là A < m \leq A + 1, hay hiểu cách khác, m chính là số nguyên nhỏ nhất lớn hơn A.
- Khi -1 \leq 0.X < 0, ta có \lfloor\frac{n}{k-1}\rfloor -1 \leq A < \lfloor\frac{n}{k-1}\rfloor, suy ra m = \lfloor\frac{n}{k-1}\rfloor.
- Khi 0 \leq 0.X < 1, ta có \lfloor\frac{n}{k-1}\rfloor \leq A < \lfloor\frac{n}{k-1}\rfloor + 1, suy ra m = \lfloor\frac{n}{k-1}\rfloor + 1.
Vậy m chỉ có thể có 1 trong 2 giá trị, hay ta chỉ cần xét m là floor hoặc ceil của phép chia n cho k-1 là đủ.