前回に引き続いてSagemathを用いて母関数から数列を求める方法の応用編です。両替問題と分割問題を扱います。
母関数を用いて両替問題を解く
両替問題は古典的な問題で母関数を作るのは簡単ですが、解くのは手計算では大変です。そこでSagemathの出番となります。以下のProject Eulerの問題を例題として取り上げます。
【Project Euler】Problem 31: コインの合計
問題の要約：1p,2p,5p,10p,20p,50p,£1(100p),£2(200p)のコインを使って£2(200p)にするのは何通りあるか求めよ
まず問題を簡略化して「1p,2p,5pコインを使って5pにするのは何通りあるか求めよ」とします。その回答は以下の式の$x^{5}$の係数となります。
(1 + x + x^{2} + x^{3} + x^{4} + x^{5}) \times
(1 + x^{2} + x^{4} ) \times
(1 + x^{5})
これをそのままSagemathで求めてみます。
n = 5
x = PowerSeriesRing(QQ, 'x', default_prec=n+1).gen()
f = (1 + x^1 + x^2 + x^3 + x^4 + x^5)*(1 + x^2 + x^4 )*(1 + x^5)
f.list()[5]
# 4
ここで1p,2p,5pのコインの母関数は既に求めたように下の表で表され、一般にnpのコインの母関数は$\Large \frac{1}{1-x^n}$なのでsagemathでは (1-x^n)^(-1)) となります。
| 母関数 | べき級数 | |
|---|---|---|
| 1p | $\Large \frac{1}{1-x^1}$ | $1 + x^1 + x^2 + x^3 + x^4 + x^5$ |
| 2p | $\Large \frac{1}{1-x^2}$ | $1 + x^2 + x^4 $ |
| 5p | $\Large \frac{1}{1-x^5}$ | $1 + x^5$ |
べき級数を母関数に変えてsagemathのプログラムにすると、以下のような簡潔なコードで表わすことが出来ました。
N = 200
coins = [1,2,5,10,20,50,100,200]
x = PowerSeriesRing(QQ, 'x', default_prec=n+1).gen()
prod([((1-x^n)^(-1)) for n in coins]).list()[N]
# 73682
母関数を用いて両替問題を解く（個数制限あり）
同じ問題の個数制限ある場合を考えます。
問題の要約： 1p,2p,5p,10p,20p,50p,£1(100p),£2(200p)のコインを使って£2(200p)にするのは何通りあるか求めよ。ただし1p,2p,5pは10個までしか使えない。
変更点以下の2点です。
- coinsをtupleのリストにして個数の上限をペアにして持たせる。上限がない場合は0。
- 母関数の生成にtruncate methodを用いて多項式の次数を制限する
Sagemathのコードは以下のようになります。2pの多項式を表示すると$x^{20}$までになっているのが分かります（なぜかオーダーが逆順に表示されますが）。
N = 200
coins = [(1,10),(2,10),(5,10),(10,0),(20,0),(50,0),(100,0),(200,0)]
x = PowerSeriesRing(QQ, 'x', default_prec=n+1).gen()
fs = [((1-x^n)^(-1)).truncate((k*n if k>0 else N)+1) for (n,k) in coins]
fs[1]
prod(fs).list()[N]
# x^20 + x^18 + x^16 + x^14 + x^12 + x^10 + x^8 + x^6 + x^4 + x^2 + 1
# 3481
母関数を用いて分割問題を解く
次にやはりProject Eulerから分割問題です。
【Project Euler】Problem 76: 整数の分割
問題の要約：100を自然数の和で何通りの組み合わせで表せるかを求めよ
分割数(Wikipedia)によると分割数は以下の母関数で表されます。
\large \prod_{k=1}^{\infty}(\frac{1}{1-x^k})
従ってそのままSagemathで実装します。
n = 100
x = PowerSeriesRing(QQ, 'x', default_prec=n+1).gen()
prod([(1-x^k)^(-1) for k in range(1,n+1)]).list()[n]
# 190569292
以上の例からわかるように、母関数が分かっている場合はSagemathに渡せばその多項式及び数列を計算してくれることが分かって頂ければ嬉しいです。
（開発環境：Cocalc）