ポアンカレ級数と SymPy で遊ぶ

この投稿は[2019年数値計算アドベントカレンダー](https://qiita.com/advent-calendar/2019/numerical_analysis)の2019年12月11日分の記事ですが、ほぼポエムです...

# ポアンカレ級数(ゆるふわな定義)
$$ V = \bigoplus_{k = 0}^\infty V_k $$
を次数付きベクトル空間 ($\deg V_k = k$) とします。このとき、$V$ の **ポアンカレ級数** $P(V; t)$ を次の形式的べき級数
$$
P(V; t) := \sum_{k = 0}^\infty (\dim V_k)\ t^k
$$
で定義します。つまり、$t$ の各べき$k$の係数が$V$の斉次成分$V_k$の次元に等しい形式的べき級数です。

## 例1. (かなり簡単な例)
$V = \mathbb{C}[x]$ (複素係数1変数多項式環)とし、$V_k, \ (k = 0, 1, 2...)$ を$k$次斉次多項式全体のなすベクトル空間としてみます。すると$\mathbb{C}[x]$ は次数付きベクトル空間と見なせますが、1変数しかないので$k$次斉次多項式は $a_k x^k, (a_k \in \mathbb{C})$ という形のものしかありません。つまり、
$$
\forall k, \ \ \dim V_k = 1
$$
となります。したがって $\mathbb{C}[x]$ のポアンカレ級数 $P(\mathbb{C}[x]; t)$ は
$$
P(\mathbb{C}[x]; t) = \sum_{k = 0}^\infty t^k = 1 + t + t^2 + t^3 + \dots
$$
と求まります。等比級数の公式を思い出すと、これは
$$
P(\mathbb{C}[x]; t) = \frac{1}{1 - t}
$$
という閉じた形に書くことができます。（形式的べき級数なので収束性は考えずにこのように書く。）

## 例2. (2変数多項式環)
今度は2変数多項式環 $V = \mathbb{C}[x, y]$ を考えてみます。$V$ の $k$次斉次成分$V_k$は上と同じように$k$次斉次多項式全体のなすベクトル空間とします。この場合は $k$ の値によって $\dim V_k$ は変わってきます。小さい方から見ていくと、

```math
\begin{align}
V_0 &= \text{Span}\langle 1 \rangle, \\
V_1 &= \text{Span}\langle x, y\rangle, \\
V_2 &= \text{Span}\langle x^2, xy, y^2\rangle, \\
V_3 &= \text{Span}\langle x^3, x^2y, xy^2, y^3\rangle, \\
&...
\end{align}
```
などとなり、$\dim V_k = k + 1$ であることがわかります。したがって
$$
P(\mathbb{C}[x, y]; t) = \sum_{k = 0}^\infty (k + 1)t^k
$$
です。これは項別微分を用いると、よい具合に変形できて、

```math
\begin{align}
P(\mathbb{C}[x, y]; t) &= \sum_{k = 0}^\infty \frac{d}{dt}t^{k + 1} \\
&= \frac{d}{dt} \left( t \sum_{k = 0} t^{k} \right) = \frac{d}{dt} \frac{t}{1 - t} = \frac{1}{(1 - t)^2}
\end{align}
```
という閉じた形に書くことができます。

## 例3. (3変数多項式環)
同様に $V = \mathbb{C}[x, y, z]$, $V_k$ は $k$次斉次多項式全体の張る部分空間とするとどうでしょうか？
$V_k$ は $x^p y^q z^r$, $p + q + r = k$ なる形の斉次式が基底となりますが、このような斉次式は

```math
\frac{(k + 2)!}{k!2!} = \dbinom{k + 2}{2}
```
種類あることがわかりますので、
$$
P(\mathbb{C}[x, y, z]; t) = \sum_{k = 0}^\infty \dbinom{k + 2}{2} t^k
$$
となります。さて、これは閉じた形に書けるでしょうか？再び項別微分を考えることでも分かりそうですが、今回は少し違う方法を考えてみます。
$$
\frac{1}{(1 - x)(1 - y)(1 - z)} = \sum_{p = 0}^\infty x^p \sum_{q = 0}^\infty y^q \sum_{r = 0}^\infty z^r
 = \sum_{p, q, r = 0}^\infty x^p y^q z^r
$$
という級数を考えてみると、右辺には $\mathbb{C}[x, y, z]$ の基底 $x^p y^q z^r$ がちょうど一回ずつ出てきています。$k$次斉次部分は $x^{p_1} y^{q_1} z^{r_1} + x^{p_2} y^{q_2} z^{r_2} + ... + x^{p_N} y^{q_N} z^{r_N}$ で $p_i + q_i + r_i = k, \ (i = 1, 2, ..., N, \ \ N = \dim(V_k))$ のようになっています。ここで $x = y = z = t$ を代入すれば $(1 + 1 + ... + 1)t^k = \dim(V_k)\ t^k$ の形にちょうどまとまるので、実は
$$
P(\mathbb{C}[x, y, z]; t) = \frac{1}{(1 - t)(1 - t)(1 - t)} = \frac{1}{(1 - t)^3}
$$
であることが分かりました。つまり
$$
\frac{1}{(1 - t)^3} = \sum_{k = 0}^\infty \dbinom{k + 2}{2} t^k \ \ \ \ \cdots \ (★)
$$
です。(ここまでで類推できるように、一般の$n$については $P(\mathbb{C}[x_1, ..., x_n]; t) = \dfrac{1}{(1 - t)^n}
$となります。)

## 数値計算（？）
さて、この(★)式が本当に成り立っていそうなことを、べき級数の最初の数項を実際に計算して確かめてみたいのですがどうすればよいでしょうか？

この記事では Python によるプログラミングで計算してみようと思います。

まず、右辺は二項係数を計算すればよいですね。Python のライブラリ [SciPy](https://www.scipy.org/) に二項係数を計算する関数 `comb` があります。
これを使って最初の10項まで計算することにすると、

```python
from scipy.special import comb
combinations = [comb(k + 2, k, exact=True) for k in range(10)]
print(combinations)
```
というプログラムを実行すればよいです。結果は

```python
[1, 3, 6, 10, 15, 21, 28, 36, 45, 55]
```
となりました。

さて、左辺についてはどうでしょうか？

$t = 0$ の周りでのテイラー展開の一般式

```math
f(t) = f(0) + f'(0)t + \frac{1}{2!}f''(0)t^2 + \frac{1}{3!}f'''(0)t^3 + ... = \sum_{k=0}^\infty \frac{1}{k!} \frac{d^kf(0)}{dt^k}t^k
```
を思い出すと、左辺の$t = 0$ における$k$階微分係数を$k!$で割ったものが係数になります。これを Python によるプログラミングで計算してみましょう。このような数式処理的な計算は [SymPy](https://www.sympy.org/en/index.html) というライブラリを使って行うことができます。同様に最初の10項まで計算することにすると**10階微分（！）**の計算が必要になり、手計算では非常に煩雑になりますが、プログラムを書くと簡単に求めることができます。
実際のプログラムは以下のようになります。(階乗の計算には `scipy`を使いました。)

```python
from sympy import *
from scipy.special import factorial
t = Symbol('t')
P = 1 / (1 - t) ** 3
coeffs = []
for k in range(10):
    dk = diff(P, t, k).subs([(t, 0)])
    coeff = dk / factorial(k, exact=True)
    coeffs.append(coeff)

print(coeffs)
```
実行結果はやはり

```python
[1, 3, 6, 10, 15, 21, 28, 36, 45, 55]
```
となり、べき級数展開の最初の10項までは左辺と右辺の係数は一致することが確かめられました。

----
# 感想
このような組合せ論的なものをコンピュータで計算したい、と日ごろ漠然と思っていたので非常に単純なものですが書いてみました。

# 参考文献
- [多項式のラプソディー (はじめよう数学)](https://www.amazon.co.jp/%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%81%AE%E3%83%A9%E3%83%97%E3%82%BD%E3%83%87%E3%82%A3%E3%83%BC-%E3%81%AF%E3%81%98%E3%82%81%E3%82%88%E3%81%86%E6%95%B0%E5%AD%A6-%E8%A5%BF%E5%B1%B1-%E4%BA%AB/dp/4535608415)