SagemathのCycle Indexを用いてピザの色分け問題を解く（その2）

前回[SagemathのCycle Indexを用いてピザの色分け問題を解く（その１）](https://qiita.com/masa0599/items/12ec8b65b7eeb3e28c19)では色の使用に制限のない場合を考えましたが。今回は制限ありの場合を考えます。母関数というものを使い、少し複雑なので一番簡単なビザを2個に切って2色を塗る場合をまず考えます。


### 【例題 2】ピザの色分け問題 (２個、2色)

これを前回と同様にの手順で考えて行きます。3.以降が変わってきます。
1. 群の置換を巡回置換表現で表す
1. 巡回指数を求める
1. 式に**色を変数**にしてを代入（母関数）
1. 母関数の多項式の係数から場合の数を求める

#### 1.群の置換を巡回置換表現で表す
巡回群C2なので以下のようになります。

| 回転置換 | 巡回置換表現 | 巡回指数の項 | 
| :--- | :--- | :--- | 
|[12] | (1)(2) | $a_1^2$| 
|[21] | (12) |$a_2$ |

#### 2. 巡回指数を求める

```math 
Z(C_2) = \frac{1}{2}(a_1^2+a_2)
```

#### 3. 式に色を変数にしてを代入

前回は色の制限がなかったので$a_i=3$と色の数を定数として代入しましたが、今回はこれを変数を使った母関数の多項式に変えます。

変数として赤の数を(r)と青の数を(b)として母関数を作るのですが、ここでは大雑把に**個数を変数の指数**、 **またはをプラス**に置き換えて多項式を作ると考えてください。

再びバーンサイドの補題のポイント **「回転して同じ形になるもの数を数えて回転の種類の数で割ったもの」** を思い出すと、例えば巡回置換$(12)$は2つの同じ色を使う必要があるので、赤が2個または青が2個。これを母関数に翻訳すると$(r^2+b^2)$になります。

色が2色場合一般に$a_i^j=(r_i+b_i)^j$になります。

| 巡回置換表現 | 巡回指数の項 | 使える色 |母関数の多項式 |
| :--- | :--- | :--- | :--- | 
| (1)(2) | $a_1^2$| $(赤が1個か青が1個)^2$| $(r+b)^2$ |
| (12) |$a_2$ | 赤が2個か青が2個|$(r^2+b^2)$ |

#### 4. 母関数の多項式の係数から場合の数を求める

それでは実際に$a_i^j=(r_i+b_i)^j$を代入してみます。

```math 
P[r,b] = \frac{1}{2}((r+b)^2+(r^2+b^2)) \\
= \frac{1}{2}(2r^2+2rb+2b^2) \\
= r^2+rb+b^2
```

これはどう解釈するかというと、組み合わせは**赤が2個、または赤と青か1個ずつ、または青が2個がそれぞれ1通りずつある**という意味になります。

この式から元の問題が「ビザを2個に切って2色を塗る場合、各色を1回ずつ使う塗り方は何通りあるか？」だとすると。答えは$rb$の係数の**1**から1通りということになります。

### Sagemathを使って解く

これをSagemathを使って解きます。巡回指数までは前回と同じです。

```python 
S = CyclicPermutationGroup(2)
S.list()
# [(), (1,2)]
P = S.cycle_index();P
# 1/2*p[1, 1] + 1/2*p[2]
```
ここに母関数の多項式を代入して展開するのですが、Sagemathではこれが以下のようにexpandを使って非常に簡単に行えます。

```python 
R.<r,b> = PolynomialRing(QQ)
pl = P.expand(2)(r,b);pl
# r^2 + r*b + b^2
```
最後に$rb$の係数を抜き出します。

```python
pl.coefficient(r*b)
# 1
```

これと同様な手順で元の問題をSagemathを使って解いていきます。

### 【例題 3】ピザの色分け問題（3色、2個ずつ）をSagemathで解く
以下のように6個に切ったピザを3色で2個ずつ塗る方法は何通りあるでしょう（回転したものは同じとみなす）

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1922159/9af1f1c6-df7c-efc1-409d-7a505efca6be.png)

巡回指数までは前回と同じです。

```python 
S = CyclicPermutationGroup(6)
P = S.cycle_index();P
# 1/6*p[1, 1, 1, 1, 1, 1] + 1/6*p[2, 2, 2] + 1/3*p[3, 3] + 1/3*p[6]
```
ここに母関数の多項式を代入して展開するのはexpandで行います。特に変数名をしていしないと自動的にx0, x1, ... が使われるようです。

```python 
p = P.expand(3);p
# x0^6 + x0^5*x1 + 3*x0^4*x1^2 + 4*x0^3*x1^3 + 3*x0^2*x1^4 + x0*x1^5 + x1^6 + x0^5*x2 + 5*x0^4*x1*x2 + 10*x0^3*x1^2*x2 + 10*x0^2*x1^3*x2 + 5*x0*x1^4*x2 + x1^5*x2 + 3*x0^4*x2^2 + 10*x0^3*x1*x2^2 + 16*x0^2*x1^2*x2^2 + 10*x0*x1^3*x2^2 + 3*x1^4*x2^2 + 4*x0^3*x2^3 + 10*x0^2*x1*x2^3 + 10*x0*x1^2*x2^3 + 4*x1^3*x2^3 + 3*x0^2*x2^4 + 5*x0*x1*x2^4 + 3*x1^2*x2^4 + x0*x2^5 + x1*x2^5 + x2^6
```
ただこれだと見づらいので.dict()を用いて指数と係数のdict形式に変換します。

```python 
pld = P.expand(m).dict();pld
# {(6, 0, 0): 1, (5, 1, 0): 1, (4, 2, 0): 3, (3, 3, 0): 4, (2, 4, 0): 3, (1, 5, 0): 1, (0, 6, 0): 1, (5, 0, 1): 1, (4, 1, 1): 5, (3, 2, 1): 10, (2, 3, 1): 10, (1, 4, 1): 5, (0, 5, 1): 1, (4, 0, 2): 3, (3, 1, 2): 10, (2, 2, 2): 16, (1, 3, 2): 10, (0, 4, 2): 3, (3, 0, 3): 4, (2, 1, 3): 10, (1, 2, 3): 10, (0, 3, 3): 4, (2, 0, 4): 3, (1, 1, 4): 5, (0, 2, 4): 3, (1, 0, 5): 1, (0, 1, 5): 1, (0, 0, 6): 1}
```
問題の答えはこのdictのkey=(2,2,2)の値となるのですが、単純にこれををkeyとしてもエラーがでてしまったので、ETupleに変換する必要があるようです。

```python
from sage.rings.polynomial.polydict import ETuple
pld[ETuple((2,2,2))]
# 16
```
無事に答えが求められてたので、これらのコードを結合して[バーンサイドの補題を用いてピザの色分け問題を解く](https://qiita.com/masa0599/items/55559687baab78254448)で作成したピザをm色でn個ずつに塗る方法の数を求める関数をSagemathの巡回指数を用いて実装した$f3(m,n)$は以下のように非常に簡単なコードになります。

```python 
from sage.rings.polynomial.polydict import ETuple
def f3(m, n):
  x = [n for i in range(m)]
  return CyclicPermutationGroup(sum(x)).cycle_index().expand(len(x)).dict()[ETuple(x)]

print(f3(3,2))
# 16
```

（開発環境：CoCalc, Sage Worksheet）