前几天又开始刷欧拉计划,遇到一道很有意思的题,把它化简之后,核心在于解这样一个组合问题:

一个袋子里有4个红球、2个蓝球,按顺序从中取出不多于2个球,问有多少种结果?(计取出顺序,且认为同色球都相同)

全部列出来就知道共有6种组合,分别是r, b, rr, rb, br, bb。推广到一般情形,即:

有 $m$ 种不同的元素,分别有 $n_1, n_2, \ldots , n_m$ 个,其中 $n=n_1+n_2+\ldots+n_m$ 。现在要从中取出不多于 $L$ 个,问有多少种可能的取法?

对于 $m=2, n_1=4, n_2=2, L=2$ 的情形,结果是6。

那对于一般的情形,又该如何计算这个值呢?

下面我将给出两种解答这个问题的思路:第一种是我当时解答这个问题的思路,第二种则是在论坛里看到的高人的解答。

解法一:动态规划算法

递推

这种解法的核心思想是递推,将一个问题转化成若干个较为简单的子问题。

记要求的取法数 $W(\overrightarrow{n}, L)$,则$W(\overrightarrow{n}=(4, 2), L=2)=6$。为了阐述方便,考虑有3个红球、2个蓝球、1个黑球的情形:

这时对第一个取出的球进行分类讨论,便可以把问题转化成三个子问题,正如下图所示:

递推关系图示

要计算原来(绿色方框中)的组合数,只需要计算下面三个(橙色方框中)的组合数即可。最终,我们便可以得到一个递推关系式:

由于每个递推关系右侧 $W(\cdot, L)$ 中的 $L$ 总是比左侧小,所以它是递推可计算的,且初值为

至此,一般情形下的递推关系式也就呼之欲出了,这里就省略不写了。

实现

还有一个问题,就是如何将上述递推关系式在计算机上实现。

一般来说,为了效率,都会将算过的函数值缓存起来,下次就可以直接使用了。因此,要使用一门带有“缓存功能”的语言来实现它。而在数据结构中,键值的缓存是用hash表来实现的。所以,若用的是C++/Java等高级编程语言,只需要用到其中的 map/HashMap 就可以了。当然 Python 也可以。实际上,只要是门编程语言,都有办法实现缓存。 (﹁ ﹁ )σ

编码

这一节是关于处理不规则键值的一个小技巧。

由于$W(\overrightarrow{n},L)$中的参数$n$是个数组/列表/向量,且维数不定,处理起来很不方便。所以我便将它进行编码,映射到一个整数。映射方式如下:

其中 $Prime[n]$ 表示第 $n$ 个质数。

首先注意到上面定义的映射是一一映射;其次它与参数的排列顺序无关;再者参数中是否包含0对结果没有影响(姑且可以认为第0个质数是1,而乘积乘1还是本身);并且编码是双向的,只要对编码值作质因数分解,就可以还原出原来的值。

这样,在实现递推函数的过程中,将编码除以一些因子,就可以给出(在递推式中)它所依赖的编码值,而这些因子可以通过质因数分解的形式得到(这句话比较绕口╮(╯-╰)╭ )。而且质因数分解得到的幂次信息,恰好代表了要求解的子问题(在递推式中)出现的次数,且自动免去了对参数列表进行排序、合并的过程(这句话自己意会去吧)。

按照上面的思路,在 Mathematica 中实现这样的过程就较为简单了。

1
2
3
4
5
6
7
8
9
rule = {3 -> 3/2, 5 -> 5/3, 7 -> 7/5, 11 -> 11/7, 13 -> 13/11};
B[L_Integer][1] = 0; (* 无元素可用 *)
B[0][k_Integer] = 0; (* 长度达到上限 *)
B[1][k_Integer] := PrimeOmega[k];(* 不同元素个数 *)
B[L_][n_] :=
B[L][n] =
Block[{fi = FactorInteger[n]},
Block[{code = n/(First /@ fi /. rule), occur = Last /@ fi},
Dot[B[L - 1] /@ code, occur] + Total@occur]]

上面的代码就描述了原来的问题,为了完整性,再给出编码函数的定义:

1
Enc[a_List]:=Times@@(a/.{i_Integer:>Prime[i]})

至此,第一种解法——动态规划 告一段落。

解法二:生成函数

概要

论坛上的高人是这么说的,“Share a interesting way”,然后贴了一段代码:

1
2
3
4
num[gred_, len_] := 
If[Length[gred] == Length[Select[gred, # >= 0 &]],
Sum[a!*Coefficient[Product[Sum[x^r/r!, {r, 0, i}], {i, gred}], x,
a], {a, 0, len}], 0];

这段代码是 Mathematica 语言中的一个函数,并且它可以直接计算出上述问题的解!

直接看想必会是一头雾水吧,因为里面涉及到一个灰常重要的概念,那就是“生成函数”。( ̄▽ ̄”)

生成函数

生成函数又称为“母函数”,如果你听说过其中任何一个名字,那想必是知道我在说什么了,如果从未听闻,也没有关系。

一个所谓的“生成函数”大概就是长的像下面这个样子:

什么呀( ⊙ o ⊙ )!这不就是 Taylor 展开麽!它就是 Taylor 展开没错,只是冠名之后就显得高大上了呀 (﹁ ﹁ )σ

为了说明它的作用,我们考虑将 10 块钱换成1块、2块、5块零钱的所有方式共有多少种。(为啥要考虑这个?大概肚子饿了想叫外卖吧,( ╯▽╰ ) 好香~~ °.°)

然后将他们都乘起来,得到

\[G(x)=(1+x+x^2+x^3+\ldots)(1+x^2+x^4+x^6+x^8+\ldots)(1+x^5+x^{10}+\ldots)\]

然后考虑 $G(x)$ 中 $x^{10}$ 项的系数,你就会惊讶地发现(也可能没发现)这恰好是将10块钱换成1块、2块、5块零钱的方法数。既然 $x^{10}$ 一定是由三个括号中每个括号里取一项乘积而成,那么这三个括号便分别对应了1块钱、2块钱、5块钱,每种货币各贡献了多少便是这个括号所提供的$x$的次数。

再进一步,有一种兑换零钱的方法,就会有一种幂次乘积的方法;反之也是如此。因此,$x^{10}$幂次前的系数就是将10块钱换成1块、2块、5块零钱的方法数。

看懂了吗?这就是“生成函数”的作用。

溯源

原问题所对应的生成函数与标准的生成函数略有不同,还是以上面“3个红球、2个蓝球、1个黑球”为例,此时的生成函数是:

有两点不同,一是生成函数变成了有限长,这是因为给的球的总数目是有限的;而是每一项都自带了一个阶乘,这是为了排除重复情形。

于是从原来的6个球中恰好取3个的方法数便是 $G(x)$ 中 $x^3$ 项的系数再乘以 $ 3! $ 。

补充解释

如下排列组合问题:

有 $m$ 种不同的元素,分别有 $n_1, n_2, \ldots , n_m$ 个,其中 $n=n_1+n_2+\ldots+n_m$ 。现在要将这 $n$ 个球排成一列,问共有多少种排列方法?

解答:分别考虑这$m$种元素的位置,即有:

若要求不超过3个的方法数,只需要求和即可,于是最终的计算公式为:

这也恰好是上面那段代码所要表达的意思!

结语

本文探讨了对于同一个组合问题的两种不同解答方式,前者是计算机的思路——动态规划,后者是数学的方法——生成函数。两种方法各有所长,且切入问题的方式完全不在同一个层面上,可最后却是殊途同归,实乃神奇之事。