呪いのかばんパズル

bagpuzzle.py

#!/usr/bin/env python3
# 呪いのカバンパズル
# https://togetter.com/li/1265539
# 以下のアルゴリズムの計算
# https://twitter.com/kyodaisuu/status/1039229105357520897
# 若干アルゴリズムを調整した

from statistics import mean, median, stdev
from scipy.optimize import fmin
import math
import sys

bag = 1000000  # 残っているカバンの数
rest = 5      # 当たっても大丈夫な回数


def put(bag, rest):
    # カバンを置く数を判定するアルゴリズム
    # @259_Momone さんの総当たり解から作成した関数
    # https://twitter.com/259_Momone/status/1039923731466866689
    # https://twitter.com/259_Momone/status/1039932917542051840
    # 平均: 36.812569 回（最小解は 36.80425）
    # 中央値: 38 回
    # 標準偏差: 5.690 回
    # 最長: 47 回
    # 最短: 5 回
    if rest == 1:
        return 1
    else:
        if bag <= 2**rest * 2.04:
            x = 2
        else:
            m = 3.478563
            n = 3.459064
            w = [1, 2.28726285, 3.61858515, 5.01286337, 6.43010308][rest-1]
            if bag <= 2**(rest * m):
                two = rest * math.log(2)
                a = (2**(rest*m/w)-2) / two**n / (m**m - 1)
                b = 2-a*two**n
                x = a * math.log(bag)**n+b
            else:
                x = bag**(1/w)
        return int(bag/x)


def E(bag, rest):
    if bag < 2:
        return 0
    p = put(bag, rest)
    hit = p/bag  # 当たる確率
    return 1 + hit * E(p, rest-1) + (1-hit) * E(bag-p, rest)


def dist(bag, rest):
    if bag == 0:
        return []
    if bag == 1:
        return [0]
    p = put(bag, rest)
    if p >= bag:
        p = bag-1
    d = []
    for i in dist(p, rest-1):
        d.append(i+1)
    for i in dist(bag-p, rest):
        d.append(i+1)
    return d


d = dist(bag, rest)
print('平均: {0} 回'.format(mean(d)))
print('中央値: {0} 回'.format(median(d)))
print('標準偏差: {0} 回'.format(stdev(d)))
print('最長: {0} 回'.format(max(d)))
print('最短: {0} 回'.format(min(d)))