演算法學習筆記：母函數詳解

引言

母函數（Generating function，生成函數）是組合數學中一種重要的方法，這裡只對最簡單的普通母函數作簡單介紹。其主要思想是，把離散序列和冪級數對應起來。

先來看一個最經典的例子：給你1克、2克、3克、4克的砝碼各一枚，問稱出1~10克的方案分別有多少種？

用母函數的方法，只需要算一個式子就好了：

$\begin{align}&(1+x)(1+x^2)(1+x^3)(1+x^4)\\=&1+x+x^2+2x^3+2x^4+2x^5+2x^6+2x^7+x^8+x^9+x^{10}\end{align}$

n次項表示稱出n克的方案，例如當n=7時有兩種方案（3+4與1+2+4）。

我們用冪級數 $\sum k_i\cdot x^i$ 表示一種狀態，它的每一項表示有 $k_i$ 種方法稱出 $i$ 克。例如 $1+x^2$ ，其實應該寫成 $1\cdot x^0+0\cdot x^1+1\cdot x^2+\cdots$ ，它的含義是，若你擁有1個2克砝碼，有1種方法稱出0克，0種方法稱出1克，1種方法稱出2克，……

當我們把兩個冪級數相乘時，我們相當於把兩種狀態的砝碼拿到了一起： $\left(\sum k_i\cdot x^i\right)\cdot\left(\sum k_j\cdot x^j\right)=\sum k_ik_jx^{i+j}$

這完全就是乘法原理，如果有 $k_i$ 種方法稱出 $i$ 克， $k_j$ 種方法稱出 $j$ 克，當然就有 $k_ik_j$ 種方法稱出 $i+j$ 克。但是 $i+j$ 可能不只有一種表示方法，我們把各種情況都加了起來，這又是加法原理。

於是，我們把看上去有點困難的組合學問題轉化為了代數問題。

我們來看例題：

（HDU1028 Ignatius and the Princess III）

“Well, it seems the first problem is too easy. I will let you know how foolish you are later.” feng5166 says.

“The second problem is, given an positive integer N, we define an equation like this:
N=a[1]+a[2]+a[3]+…+a[m];
a[i]>0,1<=m<=N;
My question is how many different equations you can find for a given N.
For example, assume N is 4, we can find:
4 = 4;
4 = 3 + 1;
4 = 2 + 2;
4 = 2 + 1 + 1;
4 = 1 + 1 + 1 + 1;
so the result is 5 when N is 4. Note that “4 = 3 + 1” and “4 = 1 + 3″ is the same in this problem. Now, you do it!”
Input
The input contains several test cases. Each test case contains a positive integer N(1<=N<=120) which is mentioned above. The input is terminated by the end of file.
Output
For each test case, you have to output a line contains an integer P which indicate the different equations you have found.

就是問一個正整數可以被用多少種方式劃分為其他正整數的和。我們可以把它轉化為：我們有無限個任意正整數重量的砝碼，有多少方式稱出n克？

我們如果有無窮多個1克砝碼，對應的母函數是什麼樣子？答案是 $1+x+x^2+x^3+x^4+\cdots$ ，也就是稱出任何自然數重量都有1種方法。相應的，如果有無窮多個2克砝碼，對應的母函數是 $1+x^2+x^4+x^6+x^8+\cdots$ 。

所以我們要算的就是 $(1+x+x^2+x^3+x^4+\cdots)(1+x^2+x^4+\cdots)\cdots$ ，這下不好辦了，無窮項的冪級數求積可不好處理。但是呢，注意到數據範圍是N<=120，那就好辦了，我們只需要考慮120次項以內的情況。

大部分題解講到這裡給出的程式，乍一看跟上面講的母函數關係不大，倒有點像DP，其實那算是優化過的，我想給出一份思路直接一些的程式碼：

#include <cstring>
#include <iostream>
using namespace std;
int pol[123][123], ans[123] = {1}, cur[123];
int main()
{
    for (int i = 1; i <= 120; ++i)
        for (int j = 0; j <= 120; j += i)
            pol[i][j] = 1; // 初始化各個多項式
    for (int i = 1; i <= 120; ++i) // ans依次乘上所有多項式
    {
        for (int j = 0; j <= 120; ++j)
            for (int k = 0; j + k <= 120; ++k) // 注意防止越界
                cur[j + k] += ans[j] * pol[i][k]; // 類似於樸素高精度乘法，還不用進位
        memcpy(ans, cur, sizeof(cur));
        memset(cur, 0, sizeof(cur));
    }
    int q;
    while (cin >> q)
        cout << ans[q] << endl;
    return 0;
}

如果不用cur，ans需要在計算過程中自己乘自己，很可能發生錯誤，所以我們用cur作為臨時數組保存中間結果。每趟多項式乘法結束後，把cur複製給ans然後清空cur。內層循環就是樸素的多項式乘法，可以用FFT、NNT優化，但這裡沒必要。

再看另一道題目：

（HDU1398 Square Coins）

People in Silverland use square coins. Not only they have square shapes but also their values are square numbers. Coins with values of all square numbers up to 289 (=17^2), i.e., 1-credit coins, 4-credit coins, 9-credit coins, …, and 289-credit coins, are available in Silverland.
There are four combinations of coins to pay ten credits:

ten 1-credit coins,
one 4-credit coin and six 1-credit coins,
two 4-credit coins and two 1-credit coins, and
one 9-credit coin and one 1-credit coin.

Your mission is to count the number of ways to pay a given amount using coins of Silverland.
Input
The input consists of lines each containing an integer meaning an amount to be paid, followed by a line containing a zero. You may assume that all the amounts are positive and less than 300.
Output
For each of the given amount, one line containing a single integer representing the number of combinations of coins should be output. No other characters should appear in the output.

就是某個國家的硬幣面額有1、4、9……元，問某個價格可以用多少種方式給出。這和剛剛那道題基本沒有什麼區別，直接計算 $(1+x+x^2+\cdots)(1+x^4+x^8+\cdots)\cdots$ ：

#include <cstring>
#include <iostream>
using namespace std;
int pol[18][305], ans[305] = {1}, cur[305];
int main()
{
    for (int i = 1; i <= 17; ++i)
        for (int j = 0; j <= 300; j += i * i)
            pol[i][j] = 1;
    for (int i = 1; i <= 17; ++i)
    {
        for (int j = 0; j <= 300; ++j)
            for (int k = 0; j + k <= 300; ++k)
                cur[j + k] += ans[j] * pol[i][k];
        memcpy(ans, cur, sizeof(cur));
        memset(cur, 0, sizeof(cur));
    }
    int q;
    while (cin >> q && q != 0)
        cout << ans[q] << endl;
    return 0;
}

砝碼有限的情形：

（洛谷P2347 砝碼稱重）

設有1g、2g、3g、5g、10g、20g的砝碼各若干枚（其總重 $\le 1000$ ），
輸入格式
輸入方式： $a_1 , a_2 ,a_3 , a_4 , a_5 ,a_6$
（表示1g砝碼有 $a_1$ 個，2g砝碼有 $a_2$ 個，…，20g砝碼有 $a_6$ 個）
輸出格式
輸出方式：Total=N
（N表示用這些砝碼能稱出的不同重量的個數，但不包括一個砝碼也不用的情況）

這道題其實完全沒必要用母函數（上面那幾道題也可以不用母函數，DP可解），但母函數也可做：

#include <cstring>
#include <iostream>
using namespace std;
int n, cnt, pol[6][1005], ans[1005] = {1}, cur[1005], step[] = {1, 2, 3, 5, 10, 20};
int main()
{
    for (int i = 0; i < 6; ++i)
    {
        cin >> n;
        for (int j = 0; n-- > -1; j += step[i]) // 循環n+1次
            pol[i][j] = 1;
    }
    for (int i = 0; i < 6; ++i)
    {
        for (int j = 0; j <= 1000; ++j)
            for (int k = 0; j + k <= 1000; ++k)
                cur[j + k] += ans[j] * pol[i][k];
        memcpy(ans, cur, sizeof(cur));
        memset(cur, 0, sizeof(cur));
    }
    for (int i = 1; i <= 1000; ++i)
        if (ans[i])
            cnt++;
    cout << "Total=" << cnt;
    return 0;
}

（這顯然可以優化，因為乘了大量的0，不過數據範圍小無所謂了）

我們也來看看其他資料給出的母函數模板是什麼意思，仍以數字劃分為例：

#include <cstring>
#include <iostream>
using namespace std;
int ans[123], cur[123];
int main()
{
    for (int i = 0; i <= 120; ++i)
        ans[i] = 1; // 第一個多項式：各項係數全部為1
    for (int i = 2; i <= 120; ++i) // 從第2個多項式開始計算
    {
        for (int j = 0; j <= 120; ++j)
            for (int k = 0; j + k <= 120; k += i) // 以當前的i為步進
                cur[j + k] += ans[j];
        memcpy(ans, cur, sizeof(cur));
        memset(cur, 0, sizeof(cur));
    }
    int q;
    while (cin >> q)
        cout << ans[q] << endl;
    return 0;
}

例如當計算 $(1+x+x^2+x^3+\cdots)(1+x^2+x^4+\cdots)$ 時，就是把它看作：

這相當於把該多項式的各次項係數分別向右移0、2、4……位，然後再按次數相加。

這種方法時間上與我的方法相差不大，但它沒有把每個多項式儲存下來，節約了空間。可以認為是利用數據特點進行的特殊優化，有一定的局限性。

母函數還有一個常用的特性：可以用等比數列求和公式把它們轉化成分式，分式相乘後再轉換回多項式。無限項的情形，可以直接假設 $0<x<1$ 然後求一個極限（例如 $1+x+x^2+\dots=\lim_{n\rightarrow\infty}\frac{1-x^n}{1-x}=\frac1{1-x}$ ）。因為這個演算法並不真的用得上 $x$ 的值，所以可以這樣假設。這樣算，有時能把問題簡化不少。例如下面這道題：

（洛谷P2000 拯救世界）

為了拯救世界，小 a 和 uim 決定召喚出 kkksc03 大神和 lzn 大神。根據古籍記載，召喚出任何一位大神，都需要使用金木水火土五種五行神石來擺一個特定的大陣。而在古籍中，記載是這樣的：
kkksc03 大神召喚方法：
金神石的塊數必須是 6 的倍數。
木神石最多用 9 塊。
水神石最多用 5 塊。
火神石的塊數必須是 4 的倍數。
土神石最多用 7 塊。
lzn 大神召喚方法:
金神石的塊數必須是 2 的倍數。
木神石最多用 1 塊。
水神石的塊數必須是 8 的倍數。
火神石的塊數必須是 10 的倍數。
土神石最多用 3 塊。
現在是公元 1999 年 12 月 31 日，小 a 和 uim 從 00:00:00 開始找，一直找到 23:00:00，終於，還是沒找到神石。不過，他們在回到家後在自家地窖里發現了一些奇怪的東西，一查古籍，哎呦媽呀，怎麼不早點來呢？這裡有一些混沌之石，可以通過敲擊而衰變成五行神石。於是，他們拚命地敲，終於敲出了n塊神石，在 23:59:59 完成了兩座大陣。然而，kkksc03 大神和 lzn 大神確實出現了，但是由於能量不夠，無法發揮神力。只有把所有用 n 塊神石可能擺出的大陣都擺出來，才能給他們充滿能量。這下小 a 和 uim 傻了眼了，趕快聯繫上了你，讓你幫忙算一下，一共有多少種大陣。

（吐槽一下，這個題面這麼長，但其實題意還表達得不是很清楚……）

要把十種召喚方法的母函數乘起來，直接做是比較困難的，但可以轉化為分式：

$1+x^6+x^{12}+\dots=\frac{1}{1-x^6}$

$1+x+x^2+\dots+x^9=\frac{1-x^{10}}{1-x}$

$1+x+x^2+\cdots+x^5=\frac{1-x^6}{1-x}$

$1+x^4+x^8+\cdots=\frac{1}{1-x^4}$

$1+x+x^2+\cdots+x^7=\frac{1-x^8}{1-x}$

$1+x^2+x^4+\cdots=\frac{1}{1-x^2}$

$1+x=\frac{1-x^2}{1-x}$

$1+x^8+x^{16}+\cdots=\frac{1}{1-x^8}$

$1+x^{10}+x^{20}+\cdots=\frac{1}{1-x^{10}}$

$1+x+x^2+x^3=\frac{1-x^4}{1-x}$

把這些分式全部乘起來，就是 $\frac1{(1-x)^5}$ 。

高數里學過冪級數展開 $(1+x)^\alpha=1+\sum_{n=1}^\infty\frac{\alpha(\alpha-1)\cdot(\alpha-n+1)}{n!}x^n$ ，用 $-x$ 替換 $x$ ，再帶進 $\alpha=-5$ ，得第 $n$ 項係數為 $\frac{(n+1)(n+2)(n+3)(n+4)}{24}$ ，注意 $(-1)^n$ 被抵消了。我們這便得到了最後的式子，只不過本題需要用高精度計算。