約瑟夫問題的兩個O(log n)解法

這個是學習編程時的一個耳熟能詳的問題了：

n個人(編號爲0,1,...,n-1)圍成一個圈子，從0號開始依次報數，每數到第m個人，這個人就得自殺，之後從下個人開始繼續報數，直到所有人都死亡爲止。問最後一個死的人的編號(其實看到別人都死了之後最後剩下的人可以選擇不自殺……)。

這個問題一般有兩種問法：

  • 給出自殺順序。不少數據結構初學書都會以這個問題爲習題考驗讀者對線性表的掌握。比較常見的解法是把所有存活的人組織成一個循環鏈表，這樣做時間複雜度是O(n*m)的。另外可以使用order statistic tree(支持查詢第k小的元素以及詢問元素的排名)優化到O(n log n)。另外有篇1983年的論文An O(n log m) Algorithm for the Josephus Problem，但可惜我沒找到下載鏈接。
  • 求出最後一個人的編號。可以套用上一種問法的解法，但另外有更加高效的解法，下文展開討論。

時間O(n)，空間O(1)

設f(n)爲初始有n人時最後一個自殺的人的編號，那麼有如下遞推式：

1
f(n) = (f(n-1) + m) mod n

以n=5, m=3爲例，一開始有這麼5個人：

1
0 1 2 3 4

第一輪報數後2號死亡，圈子裏剩下來的人的編號是這些：

1
3 4 0 1

這裏把3號寫在最前面了，因爲輪到3開始報數。如果我們有辦法知道n=4時的答案，即初始圈子爲：

1
0 1 2 3

時的答案，那麼可以把n=4的初始圈子的所有數x變換成(x+3) mod 5，得到：

1
3 4 0 1

這個恰好就是n=5時第一輪結束後的圈子狀態，也就是說我們可以根據n=4的答案推出n=5時的答案。

手工演算一下就能發現n=z時的圈子第一輪結束後(即m-1號自殺後)的圈子狀態，可以由n=z-1的初始圈子施以變換x -> (x+m) mod z得到。於是我們可以從n=1開始(此時的答案是0)，推出n=2的答案，再推出n=3，直到計算到所要求的n。

下面是C語言實現：

1
2
3
4
5
6
7
int f(int n, int m)
{
  int s = 0;
  for (int i = 2; i <= n; i++)
    s = (s + m) % i;
  return s;
}

時間O(log n)，空間O(log n)

換一個遞推思路，考慮報數過程又回到第0個人時會發生什麼。這時有floor(n/m)*m個人都已經報過數了，並且死了floor(n/m)人。之後一輪的報數會數過m個人，又會回到第0個人。

我們以n=8, m=3爲例看看會發生什麼。一開始：

1
0 1 2 3 4 5 6 7

floor(n/3)*3個人都報過數後：

1
0 1 x 3 4 x 6 7  (A)

即2號和5號死亡，接下來輪到6號開始報數。因爲還剩下6個人，我們設法做一個變換把編號都轉換到0~5：

1
2
2 3 x 4 5 x 0 1  (B)
            ___

(B)的答案等於規模縮小後的情況n'=6時最後一個人的編號s。然後根據s算出圈子狀態(A)的最後一個人的編號。

如果s在最後一個x後面(帶下劃線)，即s < n%m，那麼s2 = s-n%m+n；否則s2 = s-n%m+(s-n%m)/(m-1)；

注意如果n < m，那麼報數回到第一個人時還沒死過人。因爲子問題的規模沒有減小，所以上述方法不適用。需要用之前時間複雜度O(n)的方法遞推。下面是C語言實現：

1
2
3
4
5
6
7
int rec(int n, int m)
{
  if (n == 1) return 0;
  if (n < m) return (rec(n - 1, m) + m) % n;
  int s = rec(n - n / m, m) - n % m;
  return s < 0 ? s + n : s + s / (m-1);
}

n每次變爲n * (m-1) / m，即以一個常數因子衰減到剛好小於m，然後換用線性複雜度的遞推算法，總的時間複雜度爲O(m + log_{m/(m-1)}{n/m})，如果把m看作常數的話就是O(log n)。程序中在子問題計算後還需要做一些處理，無法尾遞歸，所以空間複雜度也等於這個。

參考：

  • Time Improvement on Josephus Problem (2002) by Fatih Gelgi
  • http://stackoverflow.com/questions/4845260/josephus-for-large-n-facebook-hacker-cup

時間O(log n)，空間O(1)

三年前我還看到過一段更巧妙的代碼，具體寫法不可考了，當時盯着式子思考了好久呢。其等價的形式長這樣：

1
2
3
4
5
6
int kth(int n, int m, int k)
{
  if (m == 1) return n-1;
  for (k = k*m+m-1; k >= n; k = k-n+(k-n)/(m-1));
  return k;
}

這個函數解決了第二種問法的一個擴展問題，即第k個(從0數起)自殺的人的編號。如果取k = n-1那麼就得到了最後一個自殺的人的編號。

這個算法的思想是第k*m+m-1次報數的人就是第k個自殺的人，追蹤他之前每次報數的時刻來確定他的編號。 考慮n=5, m=3, k=5，一共有15次報數(每報m次數死一個人，一共有k+1個人要死，所以需要(k+1)*m次報數)，每一次報數的人的編號如下：

1
2
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
0 1 2 3 4 0 1 3 4 1  3  1  3  3  3

報到2、5、8、11、14的人自殺。

設第p次(從0數起)報數的人是y，令p = m*a+b (0 <= b < m)。 經過前p次報數，一共死了floor(p/m) = a人，圈子裏還剩下n-a人。 如果y本次報數後沒有自殺，那麼他下次報數是在圈子裏其他n-a人都報數之後， 也就是第q = p+n-a = n+(m-1)*a+b次報數。 這是後繼的計算式，逆用即可計算前驅。

由y本次報數後還存活知b < m-1，故a = floor((q-n)/(m-1))， p = q-(n-a) = q-n+floor((q-n)/(m-1))。

我們要求第k個自殺的人，他是第k*m-1次報數後自殺的，用計算前驅的方式求出這個人之前一次報數是在什麼時候，不斷地重複這一過程，直到知道他是第k' (0 <= k' < n)次報數的人，那麼k'就是這個人的編號。