プロジェクトオイラー
http://projecteuler.net/index.php
Q227.
100人がテーブルの周りに並ぶ。対面の2人がそれぞれサイコロを一つずつ持つ。それぞれがサイコロを振って、1が出たら左隣、6が出たら右隣にサイコロを渡す。これが同時に行われたときに、同じ人がサイコロを2つ持ったらその時点でその人の負け。平均何回で終わるか。有効数字10桁で。
サイコロを持っている2人の距離を考える。最初は50である。1回行って、距離が変わらない確率は1/2、1変わるのは2/9、2変わるのは1/18。距離がdのときにあと何回で終わるかの期待値をE(d)とする。
18E(50) - 16E(49) - 2E(48) = 36
などとなって、50元一次連立方程式が立てられる。これを掃き出し法で解いた。
def not_zero_row(M, k):
if M[k][k] != 0:
return k
return -1def add(M, k, l):
for i in xrange(k, N + 1):
M[k][i] += M[l][i]def normalize(M, k):
a = M[k][k]
M[k][k] = 1.0
for i in xrange(k + 1, N + 1):
M[k][i] /= adef subtract(M, l, k):
a = M[l][k]
M[l][k] = 0
for i in xrange(k + 1, N + 1):
M[l][i] -= M[k][i] * aN = 50
M = [ [ 0 ] * N + [ 36 ] for n in range(N) ]
for n in range(N):
k = n + 1
if k == 1:
M[n][0] = 17.
M[n][1] = -8.
M[n][2] = -1.
elif k == 2:
M[n][0] = -8.
M[n][1] = 18.
M[n][2] = -8.
M[n][3] = -1.
elif k == N - 1:
M[n][N-4] = -1.
M[n][N-3] = -8.
M[n][N-2] = 17.
M[n][N-1] = -8.
elif k == N:
M[n][N-3] = -2.
M[n][N-2] = -16.
M[n][N-1] = 18.
else:
M[n][n-2] = -1.
M[n][n-1] = -8.
M[n][n] = 18.
M[n][n+1] = -8.
M[n][n+2] = -1.for k in xrange(N):
k0 = not_zero_row(M, k)
if k0 == -1:
exit(1)
if k0 != k:
add(M, k, k0)
normalize(M, k)
for l in xrange(k + 1, N):
subtract(M, l, k)for k in xrange(N - 1, 0, -1):
for l in xrange(k):
M[l][N] -= M[l][k] * M[k][N]print M[N-1][N]
