约瑟夫问题

约瑟夫问题的两个O(log n)解法

转自：http://maskray.me/blog/2013-08-27-josephus-problem-two-log-n-solutions
（同时添加了自己的一些心得）

这个是学习编程时的一个耳熟能详的问题了：

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

另一种说明方式， 转自 https://zh.wikipedia.org/wiki/%E7%BA%A6%E7%91%9F%E5%A4%AB%E6%96%AF%E9%97%AE%E9%A2%98

对于${\displaystyle k<n}$，可以将上述方法推广，将杀掉第k、2k、……、${\displaystyle \lfloor n/k\rfloor }$个人视为一个步骤，然后把号码改变，可得如下递推公式, , 运行时间为。

${\displaystyle f(n,k)={\begin{cases}0&{\text{if }}n=1\\(f(n-1,k)+k){\bmod {n}}&{\text{if }}1<n<k\\\left\lfloor {\frac {k((f(n',k)-n{\bmod {k}}){\bmod {n}}')}{k-1}}\right\rfloor {\text{where }}n'=n-\left\lfloor {\frac {n}{k}}\right\rfloor &{\text{if }}k\leq n\\\end{cases}}}$

这里的第三个公式，把s < n%m 与 s >= n%m，两种情况综合了起来。（公式中的k即为m）
为了理解这个公式，需要注意两个事实：

其一，下面例子中m=3，即每三个删除一个，然后从新标号得到第二行。（这里不是用约瑟夫游戏）

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

要从第二行的数得到与之对应的第一行的数，可以发现第一行的数是每三个一组，而第二行的数是每两个一组。假设第二行的一个数为v，那么v/(m-1) 可得到组号，再乘以m，即m*v/(m-1)可得到之对应的第一行的数的大概位置。然后再验算一下，发现这个结果好于预期，居然可以算出第一行的数的精确位置！

其二，约瑟夫游戏时，第二行是按递归考虑时的重新编号，

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

我们的目的是要从第二行的数得到与之对应的第一行的数，这是我们可以利用其一中的事实了，只要我们把第二行的数转为其一中的第二行的数就行了，即：

1 2	0 1 2 3 4 5 6 7 8 9 2 3 4 5 6 7 8 9 0 1

公式为 (v-2)%n, 这里n=10

-----------------------------------------------------

把上面的递归公式转为代码后，发现结果不正确，会返回一个负数。经检查，问题出在

(f(n',k)- n%k) mod n' 上。原来c++的模运算和数学意义上的还是有些差别。

from :https://stackoverflow.com/questions/7594508/modulo-operator-with-negative-values

std::cout << (-7 % 3) << std::endl;
std::cout << (7 % -3) << std::endl;

give different results:

-1
1

解释：

From ISO14882:2011(e) 5.6-4:

The binary / operator yields the quotient, and the binary % operator yields the remainder from the division of the first expression by the second. If the second operand of / or % is zero the behavior is undefined. For integral operands the / operator yields the algebraic quotient with any fractional part discarded; if the quotient a/b is representable in the type of the result, (a/b)*b + a%b is equal to a.

The rest is basic math:

(-7/3) => -2
-2 * 3 => -6
so a%b => -1

(7/-3) => -2
-2 * -3 => 6
so a%b => 1

Note that

If both operands are nonnegative then the remainder is nonnegative; if not, the sign of the remainder is implementation-defined.

就是说当f(n',k)- n%k 为负值时，(f(n',k)- n%k) mod n' 也为负。解决这个问题只需在代码中用(f(n',k)- n%k+n') mod n' 就可以了。

时间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=4，一共有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+1)*m-1次报数后自杀的，用计算前驱的方式求出这个人之前一次报数是在什么时候，不断地重复这一过程，直到知道他是第k' (0 <= k' < n)次报数的人，那么k'就是这个人的编号。

【注】对于a,b,c三个数，其中b,c确定为正，这里需要证明：
if b<c, then (a-b)/c == a/c (注意，这里不考虑小数部分)
证明：
let a-b = k*c+x,     0<=x<c
let a= k'*c+y,  0<=y<c
then
k*c+x+b=k'*c+y  then,
(k-k')c=y-b-x
here,  -b<=y-b<c-b,   -c<-x<=0, then
-2c<-b-c<y-b-x<c-b<c
but (k-k')c must be a multiple of c, so y-b-x must be 0, that means k==k'

不知道还有没有更简洁的证明方法。