\(K\)面ダイス(サイコロ)を\(N\)個投げた時に、合計がある値になる確率を求めたい
合計がある値になる場合の数を求めて\(K^N\)で割って確率にする方がやりやすそうなので、場合の数を求める方法を考える
総当り
\(N\)個のダイスについてそれぞれ\(1\)から\(K\)までループするナイーブな総当りで和を求めていくと\(O(K^N)\)ぐらいの計算量が必要になってしまい、サイコロの数が増えると現実的な時間で終わらなくなってしまう
動的計画法
漸化式を立てて動的計画法で計算すると\(O(N^2K^2)\)ぐらいで計算できる
具体的には\(n\)個サイコロを使って合計が\(c\)になる場合の数を求めたいとき\(f_n(c) = \sum_{i=1}^K f_{n-1}(c-i)\)という漸化式が成り立つ
K = 6 N = 100 dp = [[0] * N * K for i in xrange(2)] for i in xrange(K): dp[0][i] = 1 for i in xrange(1, N): for j in xrange(N * K): dp[i % 2][j] = 0 for k in xrange(max(0, j - K), j): dp[i % 2][j] += dp[(i - 1) % 2][k] for j in dp[N - 1 % 2]: print j, print
動的計画法+累積和
動的計画法の一つのサイコロに対するループで、添え字に\(1\)から\(K\)足した範囲に同じ数を足していくので、累積和に関するアルゴリズムを使えばより高速化できる
\(0\)初期化した配列を用意して、足したい範囲の開始地点に足したい数を加えて範囲の終了地点の次の添字の場所に足したい数を引いた配列を作って、各地点までの要素の累積和の配列を計算すればよい
下記記事の図がわかりやすい
計算量は\(O(N^2K)\)になる
K = 6 N = 100 dp = [[0] * N * K for i in xrange(2)] for i in xrange(K): dp[0][i] = 1 for i in xrange(1, N): dp[i % 2] = [0] * N * K for j in xrange(N * K - 1): dp[i % 2][j + 1] += dp[(i - 1) % 2][j] if j + K + 1 < N * K: dp[i % 2][j + K + 1] -= dp[(i - 1) % 2][j] cur = 0 for j in xrange(N * K): cur += dp[i % 2][j] dp[i % 2][j] = cur for j in dp[N - 1 % 2]: print j, print
フーリエ変換
確率変数の和の分布は元の分布の畳み込みになるので、フーリエ変換して係数をかけあわせて逆フーリエ変換すると目的の値が得られる(はず)
計算量は\(O(NK \log NK)\)
ただしサイコロの数を増やしたらオーバーフローしたり変な値になったりして正しい答えが得られなかった(´・ω・`)
浮動小数点演算の誤差なのかなぁ……(?)
検証してくださった方がいたので追記(2015-09-09)
@Scaled_Wurm (ΦωΦ)<このサイコロのFFTのやつ、やっぱ精度の問題で、mpmathだったらそれっぽい答えが出たです。 http://t.co/0zbfflHgTQ
— ininsanus (@ininsanus) 2015年9月9日
ソース→https://t.co/Z6qN6ebr9t
# -*- coding: utf-8 -*- import math,string,itertools,fractions,heapq,collections,re,array,bisect #http://www.slideshare.net/chokudai/fft-49066791 import cmath def pow_2_at_least(x): n = 1 while n < x: n *= 2 return n def convolution(g, N): Lg = len(g) n = pow_2_at_least(Lg * N) g = g + [0] * (n - Lg) zeta_arr = [cmath.exp(2j * cmath.pi / n) for n in xrange(1, n + 1)] gf = fft(g, n, zeta_arr) ff = [pow(gf[i], N) for i in xrange(n)] return ifft(ff, n, zeta_arr) def _fft(f, n, zeta_arr, inv): m = n / 2 if m == 1: f0 = [f[0]] f1 = [f[1]] else: f0 = [f[2 * i + 0] for i in xrange(m)] f1 = [f[2 * i + 1] for i in xrange(m)] f0 = _fft(f0, m, zeta_arr, inv) f1 = _fft(f1, m, zeta_arr, inv) sign = 1 - 2 * inv if inv: zeta = 1.0 / zeta_arr[n - 1] else: zeta = zeta_arr[n - 1] pow_zeta = 1 for i in xrange(n): f[i] = f0[i % m] + pow_zeta * f1[i % m] pow_zeta *= zeta return f def fft(f, n, zeta_arr): return _fft(f, n, zeta_arr, False) def ifft(f, n, zeta_arr): return [i / n for i in _fft(f, n, zeta_arr, True)] K = 6 N = 100 print '\n'.join([str(int(round(i.real))) for i in convolution([0] + [1.0] * K, N)][1:N * K + 1])
参考URL
いろいろな確率とかが載ってる
\(K\)面ダイスを\(N\)個投げた時に合計がある値\(C\)になる確率を求める以下の式の導出が載っている
$$P(N, K, C)=\frac{1}{K^N}\sum_{i=0}^{\lfloor \frac{C - N}{K} \rfloor}(-1)^i \binom Ni \binom{C-si-1}{N-1}$$