205. 骰子游戏(Dice Game)

彼得有九个四面(金字塔形)骰子，每个面的编号分别为1、2、3、4。科林有六个六面的(立方体)骰子，每个面的编号分别为1、2、3、4、5、6。彼得和科林掷骰子并比较各个骰子编号的总和，总和高的人胜利； 如果总和相等，则结果为平局。请问金字塔彼得击败立方体柯林的概率是多少(以0.abcdefg的形式将答案四舍五入到小数点后七位)？

分析：这道题要求彼得赢的概率，即其抛出的骰子面值之和大于科林抛出的骰子面值之和，假设彼得每次抛出骰子的面值之和为\(ps_i\)，科林每次抛出的骰子面值之和\(cs_i\)，则彼得获胜的概率\(P\)为：

彼得和科林抛出特定面值和的概率应为其抛出这个面值的可能性总数除以所有可能性总数，因为彼得的总共有九个四面骰子，则其可能性总数为\(4^9\)种；同理科林有六个六面骰子，所以总的可能性为\(6^6\)种。设彼得抛出特定面值和的可能性数为\(pn_i\)，科林为\(cn_i\)，则有：

因此，为计算出彼得获胜的概率，我们需要计算出彼得与科林抛出特定面值和可能性总数。这个问题可以一般化为对于\(m\)个\(n\)面的骰子，其面值之和为\(s\)的可能事件总数是多少，我们可以将其设为一个函数，也即求\(N(m,n,s)\)的值。为了求这个值，至少有两种方法：

一、暴力枚举法

考虑到这个问题的规模不大，我们可以直接用枚举的方式得到问题的解，基本思路是：使用\(itertools\)模块中\(product\)函数，枚举所有骰子可能出现的面值，分别计算不同情况下骰子的面值之和，最后使用\(colletions\)模块中的\(Counter\)类统计不同面值出现的次数，结果如下：

从上表我们可以看出，彼得能够抛出的九个四面体的面值之和最小为9，最大为36，可以抛出的次数都是一次；而科体能够抛出的六个六面体的面值之和最小为6，最大为36，可以抛出的次数也都是一次。假设现在彼得抛出的面值之和为9，为确保其比科林抛出的大，那么科林抛出的面值之和只能为6、7、8，则此时彼得获胜的概率为：

二、生成函数法

上面这种方法在问题规模较小时可以比较好的应对，也相对更容易理解。但是如果问题规模进一步增加，比如色子的数量达到一百，那么这种暴力枚举的方法就力不从心了。一种更为高端的方法的使用生成函数的方法求出函数\(N(m,n,s)\)的解析式，从而直接计算出相应的值。具体的方法大家可以参考这篇论文，我这里只给出一个结果：

比如对于九个四面体色子，它们面值之和等于15的可能性总数可以如下求解：

和上面表格显示的结果一致。

三、结论

对于这道题的数据规模，以上两种方法都可以在较快时间得到答案，第二种方法约比第一种方法快三十倍。但对于数据规模更大的问题，第二种方法显然效率会更高一些。两种方法我都写了相应的代码，具体如下：

# approach 1, time cost = 101 ms ± 484 µs

from itertools import product
from collections import Counter

def dice_values_sum_ways(sides,number):
    arrangement = product(range(1,sides+1),repeat=number)
    sum_of_values = [sum(x) for x in arrangement]
    c = Counter(sum_of_values)
    return dict(c)

def main():
    beats = 0
    total = (4**9) * (6**6)
    pw = dice_values_sum_ways(4,9)
    cw = dice_values_sum_ways(6,6)
    for k1,v1 in pw.items():
        for k2,v2 in cw.items():
            if k1 > k2:
                beats += (v1*v2)
    return round(beats/total,7)

# appraoch 2, time cost = 3.88 ms ± 53.3 µs

from scipy.special import comb

def numbers_of_ways(m,n,p):
    up = int((p-m)/n)+1
    seq = [(-1)**j * comb(m,j) * comb(p-n*j-1,m-1) for j in range(up)]
    return sum(seq)

def prob():
    beats = 0
    total = (4**9) * (6**6)
    pw = {i:numbers_of_ways(9,4,i) for i in range(9,37)}
    cw = {i:numbers_of_ways(6,6,i) for i in range(6,37)}
    for k1,v1 in pw.items():
        for k2,v2 in cw.items():
            if k1 > k2:
                beats += (v1*v2)
    return round(beats/total,7)