欧拉计划：一个组合问题

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

一个袋子里有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((3,2,1),3)=W((2,2,1),2)+W((3,1,1),2)+W((3,2),2)+3$$

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

$$W(\overrightarrow{n},1)=\text{dimension of }n$$

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

实现

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

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

编码

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

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

$$Enc[(a_1,a_2,\ldots,a_s)]=Prime[a_1] \times Prime[a_2] \times\ldots\times Prime[a_s]$$

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

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

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

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

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]]

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

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

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

解法二：生成函数

概要

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

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 语言中的一个函数，并且它可以直接计算出上述问题的解！

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

生成函数

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

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

$$g(x)=\frac{1}{1-x}=1+x+x^2+x^3+\ldots$$

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

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

记

$$\begin{cases} g_1 (x)&=\frac{1}{1-x}=1+x+x^2+x^3+\ldots \\ g_2 (x)&=\frac{1}{1-x^2}=1+x^2+x^4+x^6+x^8+x^{10}+\ldots \\ g_5 (x)&=\frac{1}{1-x^5}=1+x^5+x^{10}+\ldots \end{cases}$$

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

\[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个黑球”为例，此时的生成函数是：

$$\begin{cases} g_{red}(x)&=1+x+\frac{x^2}{2!}+\frac{x^3}{3!} \\ g_{blue}(x)&=1+x+\frac{x^2}{2!}\\ g_{black}(x)&=1+x \end{cases}$$

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

于是从原来的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$种元素的位置，即有：

$$M=C_{n}^{n_1} C_{n-n_1}^{n_2} \ldots C_{n_m}^{n_m}=\frac{n!}{n_1! n_2! \ldots n_m!}$$

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

$$W(n, L)=\sum _{a=1}^L a! \times (\text{the coefficient of }x^a\text{ in }\prod _{i\in n} \sum _{k=0}^i \frac{x^k}{k!})$$

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

结语

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