囚人問題

きり@kiri8128 さん出題の問題
https://twitter.com/kiri8128/status/1297582612668538882

午前零時がどの日に属しているかは曖昧なので、スイッチをオン/オフにした日の夜に、それに応じた照明の点灯が行われると考えることにする。

0° 基本的なアイデア

囚人の集合 P が与えられていて、各囚人は P に自分が属しているどうかを知っているとする。
また、P に属さない囚人の人数は M 人以下であることが判明しているとする。
このとき、P=\emptyset かどうかを以下の方法で囚人全員が知ることができる。

初日に P に属する囚人がスイッチを入れる。
2日目から M 日目まで、前日にスイッチを入れた、または自分のいた部屋の照明が点灯した囚人はスイッチを入れる。

スイッチを入れたことがある囚人と入れたことのない囚人が共存していたとすると、スイッチを入れたことのない囚人の人数は1日で少なくとも一人減るので、P\neq\emptyset なら全ての部屋の照明が点灯し、P=\emptyset なら M 日目の夜にどの部屋の照明も点灯しない。

1° 囚人数の上限を調べる

K=1,2,3,\cdotsK を増やしながら以下を行う。

初日は「あなた」一人がスイッチを入れる。
2日目から K+1 日目まで、前日にスイッチを入れた、または自分のいた部屋の照明が点灯した囚人はスイッチを入れる。
K+1 日目の夜に自分の部屋の照明が点灯した囚人の人数は高々 2^K 人。
K+1 日目の夜に自分の部屋の照明が点灯しなかった囚人の集合を PM=2^K として、0° の手続きによりP=\emptyset であるかどうかを確かめる。

以上を繰り返して P=\emptyset となった時点で終了。この時点で囚人の人数が 2^K 人以下であることを全員が知ることができる。
以降の議論ではこの K が判明していることを仮定する。

2° 囚人の人数を求める

囚人全体を排他的な m 個の集合に分割する。
初期状態は、m=2 で、「あなた」ひとりから成る集合を G_1、それ以外の囚人の集合を G_2 として、 |G_1| = 1, |G_2| = x_1 とする。x_1 は未知のパラメーター。
必要ならこの分割をより細かくしながら、未知のパラメーターを消去していく。未知のパラメーターが全て消去できれば、囚人の人数が確定する。

以下のステップを繰り返す。

分割された囚人の集合を  G_1, G_2, \ldots, G_m とする。
 a_i, c_{ij} を定数、未知のパラメーターを  x_1, x_2, \ldots, x_n とする。
 |G_i| = a_i + \sum_{j=1}^n c_{ij} x_j\ \ \ (i=1,2,\ldots,m)
であることが分かっているとする。
囚人の人数が確定していないなら  c_{ij}\neq 0 なる係数が存在するはずなので、そのうちのひとつを選ぶ(どう選ぶかのアルゴリズムは最初のメールで囚人全員が共有しておく)。
 c_{ij} c_{kj} > 0 なる全ての k(k=1,2,\ldots,m) について  H=\bigcup G_k とする。
H に属する囚人がスイッチを入れて、その夜に部屋の照明が点灯した囚人の集合を S とし、S に属さない囚人の集合を {\bar S} とする。
G_l について P=S\cap G_l, M=2^K として 0° の手続きにより S\cap G_l =\emptyset であるかどうかを確かめる。{\bar S}\cap G_l=\emptyset についても同様に確かめる。
 S \cap G_l\neq \emptyset かつ  {\bar S} \cap G_l\neq\emptyset なる各 G_lG'_l=S\cap G_lG''_l={\bar S} \cap G_l のふたつの集合に分割し、新しいパラメーター y_l を導入して  |G'_l|=y_l, |G''_l| = |G_l| - y_l とする。

|S||S| = \sum_{l=1}^m |S\cap G_l| で、和の中の各項は S\cap G_l=\emptyset なら当然 |{\bar S}\cap G_l|=0{\bar S} \cap G_l=\emptyset なら |S\cap G_l|=|G_l|、どちらでもなければ |S\cap G_l|=|G'_l|=y_l で計算できる。
スイッチを入れた囚人の人数と照明が点灯した部屋の囚人の人数は同じであることから |H|=|S| が成立し、パラメーターについての方程式が一本定まる。この方程式によってパラメーター  x_j を消去する(H の選び方により x_j は必ず消去できるのだが、その理由は下に書く)。

集合が分割されるたびにパラメーターが1個導入され、分割された集合の数は囚人の人数 N 以下なので、導入されるパラメーターは高々 N-2 個(初期状態で集合が2個だから)。初期状態でパラメーターが1個あるので消去されるべきパラメーターは高々 N-1 個。1ステップごとにひとつのパラメーターが消去されるので、この操作は N-1 以下のステップ数で終了する。

3° 具体例

(1) 最初、「あなた」ひとりから成る集合を G_1、それ以外の囚人の集合を G_2 として、 |G_1| = 1, |G_2| = x_1
H=G_2 として S を定めて、0° の手続きにより
 {\bar S}\cap G_1=0, S\cap G_2\neq\emptyset,{\bar S}\cap G_2\neq\emptyset
が分かったとする。
G_2G'_2 = S\cap G_2G''_2={\bar S} \cap G_2 に分割して、パラメーター y_2 を導入して、|G'_2|=y_2, |G''_2|=x_1-y_2
|H|=|G_2|=x_1
|S|=|S\cap G_1| + |S\cap G_2| = |G_1| + |G'_2|=1+y_2
({\bar S}\cap G_1=\emptyset より S\cap G_1 = G_1 であることを使った。)
|H|=|S| より x_1=1+y_2
これにより  |G''_2| から x_1 を消去して |G''_2| = 1

(2) 集合とパラメーターの名前を付け直して、|G_1|= |G_2|=1, |G_3|=x_1
H=G_3 として S を定めて、
 {\bar S}\cap G_1={\bar S}\cap G_2=\emptyset,S\cap G_3\neq\emptyset, {\bar S}\cap G_3\neq\emptyset
が分かったとする。
G_3G'_3 = S\cap G_3G''_3={\bar S} \cap G_3 に分割して、パラメーター y_3 を導入して、|G'_3|=y_3, |G''_3|=x_1-y_3
|H|=|S| より x_1=2+y_3
これにより  |G''_3| から x_1 を消去して |G''_3| = 2

(3) 集合とパラメーターの名前を付け直して、|G_1|=|G_2|=1, |G_3|=2, |G_4|=x_1
H=G_4 として S を定めて、
 {\bar S}\cap G_1=S\cap G_2={\bar S}\cap G_3=0, S\cap G_4\neq\emptyset, {\bar S}\cap G_4\neq\emptyset
が分かったとする。
G_4G'_4 = S\cap G_4G''_4={\bar S} \cap G_4 に分割して、パラメーター y_4 を導入して、|G'_4|=y_4, |G''_4|=x_1-y_4
|H|=|S| より x_1=3+y_4
これにより  |G''_4| から x_1 を消去して |G''_4| = 3

(4) 集合とパラメーターの名前を付け直して、|G_1|=|G_2|=1, |G_3|=2, |G_4|=3, |G_5|=x_1
H=G_5 として S を定めて、
 S\cap G_1=S\cap G_2=S\cap G_3=\emptyset,
S\cap G_4\neq\emptyset, {\bar S}\cap G_4\neq\emptyset, S\cap G_5\neq\emptyset, {\bar S}\cap G_5\neq\emptyset
が分かったとする。
G_4G'_4 = S\cap G_4G''_4={\bar S} \cap G_4 に分割して、パラメーター y_4 を導入して、|G'_4|=y_4, |G''_4|=3-y_4
G_5G'_5 = S\cap G_5G''_5={\bar S} \cap G_5 に分割して、パラメーター y_5 を導入して、|G'_5|=y_5, |G''_5|=x_1-y_5
|H|=|S| より x_1=y_4+y_5
これにより  |G''_5| から x_1 を消去して |G''_5| = y_4

(5) 集合とパラメーターの名前を付け直して、|G_1|=|G_2|=1, |G_3|=2, |G_4|=|G_5|=x_1, |G_6|=3-x_1, |G_7|=x_2
H=G_4 \cup G_5 として S を定めて、
 {\bar S}\cap G_1=S\cap G_2=S\cap G_3=S\cap G_4=S\cap G_5=S\cap G_6={\bar S}\cap G_7=\emptyset
が分かったとする。
|H|=|S| より 2x_1=1+x_2
これにより  |G_4|, |G_5|, |G_6| から x_1 を消去して |G_4|=|G_5|= \frac{1+x_2}{2}, |G_6|=\frac{5-x_2}{2}

(6) パラメーターの名前を付け直して、|G_1|=|G_2|=1, |G_3|=2, |G_4|=|G_5|=\frac{1+x_1}{2}, |G_6|=\frac{5-x_1}{2}, |G_7|=x_1
H=G_7 として S を定めて、
 {\bar S}\cap G_1=S\cap G_2=S\cap G_3=S\cap G_4=S\cap G_5=S\cap G_6=S\cap G_7=\emptyset
が分かったとする。
|H|=|S| より x_1=1
これにより  |G_4|, |G_5|, |G_6|, |G_7| から x_1 を消去して |G_4|=|G_5|=|G_7|=1, |G_6|=2

以上でパラメーターが全て消去され、囚人の人数は
 N= \sum_{i=1}^7 |G_i| = 9
と確定する。

パラメーターが消去できること

2° でパラメーター x_j が必ず消去できると書いたが、上の (5) のステップを例に説明する。
(5) では |G_1|=|G_2|=1, |G_3|=2, |G_4|=|G_5|=x_1, |G_6|=3-x_1, |G_7|=x_2 となっていたわけだが、例えば H=G_4 とかとしてしまうと、看守に意地悪されてパラメーターが消去できないことがある。
H から S を定めて {\bar S}\cap G_5 = \emptyset
S\cap G_1 = S\cap G_2 = S\cap G_3 = S\cap G_4 = S\cap G_6 = S\cap G_7=\emptyset となったとしよう(つまり G_5 の部屋だけ照明が点灯したとする)。
このとき、|H|=|S|x_1=x_1 となって、この式からは x_1 が消去できなくなってしまう(つまり有効な情報が得られていない)。

2° に書いたような H の選び方をすれば、こういう不都合なことが起こらないことを説明する。
2° で書いたのは要するに、適当なパラメーター x_j を選んで、|G_k| のパラメーター x_j の係数が例えば正になる G_k だけを全て拾ってきて H=\bigcup G_k とせよということだった。
(5) の場合、x_1 の係数が正のものを拾うことにして、 該当するのは G_4, G_5 だけだから H=G_4 \cup G_5 としたのだった。
このとき、|H|=2 x_1 だから |H|x_1 の係数は 2。
|S|=\sum_{l=1}^7 |S\cap G_l| で、この和のうち x_1 の係数が正になり得るのは |S\cap G_4| + |S\cap G_5| の部分だけ。
|S\cap G_4| が取り得る値は、囚人の配置によるが x_1, 0, y_4 のうちのどれか。|S\cap G_5| についても同様。
つまり、S\cap G_4=G_4, S\cap G_5=G_5 のときだけ |S|x_1 の係数が 2 になる可能性があって、それ以外は 2 未満にしかならない。2 未満なら不都合は起こらないので、S\cap G_4=G_4, S\cap G_5=G_5 となった場合だけが問題になる。
このとき、|S|\geq |S\cap G_4| + |S\cap G_5| = |G_4| + |G_5|
H=G_4 \cup G_5, |H|=|S| と併せて考えれば結局 H=S
H=S というのはスイッチを入れた囚人の部屋の照明が全て点灯するということだから、H の囚人が円環状に配置されていなければならない。しかし「あなた」は H には属さないのでこういうことはあり得ない(今の例に限らず、2° の方法に従う限り、一般の場合も「あなた」は H に属さない)。
結局、H=G_4 \cup G_5 なら不都合なことは起こらない。

今は (5) を例に取って考えたが、一般の場合も(たぶん)明らか。

「メダルツアー」問題

@riverplus さん出題の問題。
問題: http://riverplus.hatenablog.com/entry/2018/02/25/143351
解説: http://riverplus.hatenablog.com/entry/2018/03/19/220532

下から k 番目 (0≦k<h) のメダルを M_k、C(n,k) を2項係数として、
スタートから M_k の下の点までの経路の数は C(m + k - 1, k)、
M_k の上の点からゴールまでの経路の数は C(w + h - m - k, h - k - 1) なので、
M_k を拾う経路の数は C(m + k - 1) C(w + h - m - k, h - k - 1)。
M_k を拾う確率は、これを経路の総数 C(w + h, w) で割ればよい。
拾うメダルの個数の期待値 F(w,h,m) は、k=0 から h-1 まで M_k を拾う確率を足し上げればよい。
\Large F(w,h,m) = \frac{1}{C(w+h,w)} \sum_{k=0}^{h-1} C(m+k-1,k)C(w+h-m-k,h-k-1)

Ruby で書くとこうなる。

@fac = [1]
def fac(n)
  @fac[n] ||= n * fac(n - 1)
end

def C(n, k)
  fac(n) / fac(k) / fac(n - k)
end

def F(w, h, m)
  (0...h).map {|k| C(m + k - 1, k) * C(w + h - m - k, h - k - 1)}.inject(:+) / C(w + h, w).to_f
end

DATA.each do |s|
  w, h, m = s.split.map(&:to_i)
  puts "F(%d, %d, %d) = %.7f" % [w, h, m, F(w, h, m)]
end

__END__
3 2 2
1 1 1
2 1 3
4 3 2
10 10 3

以上で解決と思っていたら、解説記事にあるとおり、ずっと簡単な解法があった。
上に一歩歩くのを 'U'、右に一歩歩くのを 'R' と書くと、U を h 個、R を w 個並べる方法の総数が経路の総数 C(w + h, w) になる。
これは、U を h 個並べておいて、その(両端を含めた) h+1 個の隙間に w 個の R を詰め込む方法の総数に等しい。つまり h+1 部屋に w 人を押し込む方法の数と同じ。
F(m, w, h) は、m 番目の部屋の人数の期待値に等しく、h+1 個の部屋は全部平等なので、
\Large F(m, w, h) = \frac{w}{h+1}
となる。

辞書順に数を並べる

@Nabetani さんの問題
http://nabetani.sakura.ne.jp/hena/orde22numord/
を解いてみた。
m 以上 n 以下の自然数を b 進表記して辞書順に並べたとき、x 番目の数を求めよという問題。
他の方の実装が
http://qiita.com/Nabetani/items/99e38a39165e1905b415
にある。

だらだらやってたら3時間くらいかかった。
言語は Ruby

# u を kn 桁(= n の桁数)の数として、順序が u の順序以下のものの数
def f(u, m, n, b, km, kn)
  v = b * u - b ** kn
  [u / b ** (kn - km) - m + 1, 0].max - [u - n, 0].max +
  (km...kn).inject(0) {|s| s + (v /= b) + 1}
end

def solve(m, n, b, x)
  km, kn = [m, n].map {|j| j.to_s(b).size}
  u = (b ** (kn - 1)...b ** kn).bsearch {|u| f(u, m, n, b, km, kn) >= x}
  # この時点で、u は目的の数の末尾に0個以上の '0' が付いている
  u / b ** (f(u, m, n, b, km, kn) - x + (u > n ? 1 : 0)) # 末尾の 0 を削る
end

DATA.each do |s|
  j, t, e = s.split
  puts "%2s %-30s %d" % [j, t, ans = solve(*t.split(',').map(&:to_i))]
  fail unless ans == e.to_i
end

__END__
0   2,15,3,8                        14
1   6,8,8,1                         8
2   6,6,5,1                         6
3   4,11,5,2                        6
4   9,17,8,2                        10
5   73,82,18,5                      77
6   70,149,2,48                     95
7   82,119,15,26                    107
8   40,127,26,47                    86
9   851,950,31,89                   939
10  660,807,34,143                  802
11  962,1227,34,186                 1075
12  544,1258,25,245                 869
13  1208,7380,2,803                 4630
14  2,100000000,10,1                10
15  1642,6626,3,4354                2029
16  8623,15873,2,4188               12810
17  7013,16409,17,3919              10931
18  7222,13243,19,3623              10844
19  4421,87412,5,62080              60697
20  51812,61593,22,2957             54768
21  44260,67742,17,15128            59387
22  90929,696974,36,84356           175284
23  1,100000000,2,50000000          92108858
24  101262,924931,7,341846          358105
25  211114,675661,33,312769         523882
26  1,100000000,36,50000000         11855060
27  1,100000000,10,100000000        99999999
28  826566,1444644,36,186226        1012791
29  6687091,6985458,4,251580        6938670
30  4219082,6179401,12,538113       4757194
31  7781931,8634872,23,217457       7999387
32  818218,64845924,12,57752258     30039084
33  4987258,11760092,32,6741567     11728824
34  7075964,28431054,22,11548591    18624554
35  26065964,66252692,18,29196596   63208819
36  66097362,104618756,16,32740764  98838125

数独ソルバー

Ruby だと1行のメソッドで書ける気がしたので書いてみた。
効率は悪いけど、大体数秒から数分で解ける。
問題は1行の文字列で与える。'0' が空白。
解答は可能な解を配列に並べたもの(解が複数あるときは全部求める)。

例題は Dr. Arto Inkala 作の世界一難しいらしいという問題を借りてきた。この問題だと数秒で解ける。

def solve(s, k = 0)
  (k = (k..80).find {|k| s[k] == '0'}) ? (('1'..'9').to_a - (0..8).flat_map {|i|
    [s[k/9*9 + i], s[k%9 + i*9], s[k%9/3*3 + k/27*27 + i + i/3*6]]
  }).flat_map {|c| solve(s[0...k] + c + s[k + 1..80], k + 1)} : [s]
end

p solve('005300000800000020070010500400005300010070006003200080060500009004000030000009700')

Berlekamp-Massey のアルゴリズム

線形の漸化式(整数係数で最高次の係数は 1 とする)に従う数列が与えられたとき、もとの漸化式の特性方程式を見つける方法があるというのを某所で知った。
例えば、フィボナッチ数列の何項か(例えば 1, 2, 3, 5, 8)が与えられたとき、特性方程式
x^2 - x - 1=0
を求めるとかそういうこと。
これには Berlekamp-Massey のアルゴリズムというのを使えばいいらしい。

Berlekamp-Massey のアルゴリズムのコードは英語版 Wikipedia にあるので、適当に Ruby に直してみた。
Berlekamp-Massey のアルゴリズムは任意の体で使えるらしいのだが、有理数体で計算すると途中で巨大な分子・分母の分数が出てきて計算が進まなくなるので、いくつかの素数 p について GF(p) で求めた多項式から、中国剰余定理を使って、有理数体での特性方程式を求めるようにしている。
この部分は自分で考えているので変なことをしているかもしれない。

以下のコードは入力する数列が短かすぎると無限ループになってしまうようなので、その場合は数列を長くして再計算する必要がある。
数列の長さが特性方程式の次数(この次数があらかじめ分かっていないこともある)の2倍以上ならうまくいくようだが、必ずうまくいくかどうかは知らない。

require 'prime'

# 多項式は係数を並べた配列で表現する
# 例えば x^2 + 2x + 3 は [1, 2, 3] と表す

# Berlekamp-Massey のアルゴリズム
# q を素数として GF(q) で計算し、特性方程式 c(x)=0 の c を返す
def berlekamp_massey(s, q)
  b, c = [1], [1] + [0] * (s.size - 1)
  l, m, a = 0, -1, 1
  s.size.times do |n|
    d = (0..l).inject(0) {|sum, i| (sum + c[i] * s[n - i]) % q}
    next if d == 0
    t = c[0..l]
    (0...[s.size - n + m, b.size].min).each do |j|
      c[n - m + j] = (c[n - m + j] - d * a * b[j]) % q
    end
    b, l, m, a = t, n + 1 - l, n, mod_inv(d, q) if 2 * l <= n
  end
  c[0..l]
end

def euclid(a, b)
  return [0, 1] if a == 0
  q, r = b.divmod(a)
  x, y = euclid(r, a)
  [y - q * x, x]
end

# x^(-1) (mod n)
def mod_inv(x, n)
  euclid(x, n)[0]
end

# x % n1 = r1, x % n2 = r2, |x| <= n1 * n2 / 2 となる x
def chinese(n1, r1, n2, r2)
  x = (n1 * (r2 - r1) * mod_inv(n1, n2) + r1) % (n = n1 * n2)
  2 * x > n ? x - n : x
end

# f を多項式として f=0 が数列 s を生成する漸化式の特性方程式となっているか
def test(f, s)
  (0..s.size - f.size).all? do |i|
    f.each_with_index.inject(0) {|sum, (fj, j)|
      sum + fj * s[f.size + i - j - 1]
    } == 0
  end
end

# 数列 s を生成する漸化式の特性方程式 f=0 の f を返す
def polynomial(s)
  f, n = [], 1
  Prime.each do |q|
    c = berlekamp_massey(s, q)
    if c.size != f.size then
      f, n = c, q if c.size > f.size
      next
    end
    f = (0...f.size).map {|i| chinese(n, f[i], q, c[i])}
    return f if test(f, s)
    n *= q
  end
end


# 使用例 1
# フィボナッチ数列の特性方程式を求める
p polynomial([1, 2, 3, 5, 8])


# 使用例 2
# フィボナッチ数とトリボナッチ数の積の数列の特性方程式を求める

# フィボナッチ数を F_i, トリボナッチ数を T_i として、s_i = F_i * T_i とする
# n は項数
def init_s(n)
  f, t = [0, 1], [0, 0, 1]
  n.times.each_with_object([]) do |i, s|
    s << f[-1] * t[-1]
    f << f[-1] + f[-2]
    t << t[-1] + t[-2] + t[-3]
  end
end

p polynomial(init_s(24))