SageMathでの連立方程式の解き方(多項式環含む)

# I. はじめに

[corCTF 2022](https://2022.cor.team/)のrlfsrという問題がありまして，連立方程式を解く際にどうアプローチするか(やコーディング方法)に色々と悩んだため備忘録として残します．

# II. 実数での解法

## II-a. 行列形式で表記できる場合

$$
\begin{cases}
x + 2y + 3z = 0 \\\\
3x + 2y + z + 4 = 0 \\\\
x + y + 2z + 1 = 0
\end{cases}
$$

書き換えると

$$
\begin{pmatrix}
1 & 2 & 3 \\\\
3 & 2 & 1 \\\\
1 & 1 & 2
\end{pmatrix}
\begin{pmatrix}
x \\\\
y \\\\
z
\end{pmatrix}=
\begin{pmatrix}
0 \\\\
-4 \\\\
-1
\end{pmatrix}
$$

解は $(x, y, z) = (-2, 1, 0)$です．

```python
A = Matrix([[1, 2, 3],
            [3, 2, 1],
            [1, 1, 2]])
v = vector([0, -4, -1])

print(A.solve_right(v))
```


## II-b. 行列形式では表記できない場合(変数が少ない場合)

これは非推奨．II-c.の方がよいです．

$$
\begin{cases}
xy + z = 0 \\\\
x + yz = 0 \\\\
x + y + z + 1 = 0
\end{cases}
$$

解は $(x, y, z) = (0, -1, 0)$です．

右辺の指定がない場合はデフォルトで $ = 0$ として扱われるので省略しても問題ないです．

```python
var('x y z')

eq1 = x * y + z  # "eq1 = x * y + z == 0" is also OK.
eq2 = x + y * z
eq3 = x + y + z + 1

print(solve([eq1, eq2, eq3], x, y, z))
```

## II-c. 行列形式では表記できない場合(変数がそこそこ多い場合)

ほぼII-b. と同じですが，例えば30変数を扱いたい場合に，II-b. の1行目で `var(x_1, ..., x_30)`とはさすがに書いていられないです…．

なので多項式環を使います．

例としてはII-b. の

$$
\begin{cases}
xy + z = 0 \\\\
x + yz = 0 \\\\
x + y + z + 1 = 0
\end{cases}
$$

を流用します．

```python
PR = PolynomialRing(QQ, 3, "x")
xvars = PR.gens()

# print(xvars) => (x0, x1, x2)

eq1 = xvars[0] * xvars[1] + xvars[2]
eq2 = xvars[0] + xvars[1] * xvars[2]
eq3 = xvars[0] + xvars[1] + xvars[2] + 1

J = Ideal([eq1, eq2, eq3])
print(J.variety())
```

「II. 実数の場合」と言っておきながら1行目で`QQ`として有理数に限定していますが，このほうが高速に動く気がします．
(無理数を扱うCTFの問題がそもそも少ないですし)

もしどうしても実数で扱いたい場合は`QQ`ではなく`RR`としましょう．


## II-d. 行列形式では表記できない場合(変数がとても多い場合．解決策未発見)

以下の67変数の式を考えてみます．

$$
\begin{cases}
x_0 + 1 = 0 \\\\
x_1 + 1 = 0 \\\\
\vdots \\\\
x_{66} + 1 = 0
\end{cases}$$

当然ながら解は $x_0 = x_1 = \cdots = x_{66} = -1$ です．

これをII-c. 通りに実装すると以下のようになりますがエラーを吐きます．

```python
PR = PolynomialRing(QQ, 67, "x")
xvars = PR.gens()

eqs = [_x + 1 for _x in xvars]

J = Ideal(eqs)
print(J.variety())
```

長いのでエラーコードは貼りませんが，一番最後に

```
RuntimeError: error in Singular function call 'groebner':
long int overflow in hilb 4
```

とあるのでキャパオーバーでしょう．

このようなシンプルな場合でさえ67変数で動かないので(66変数だと動作確認済)，条件が複雑だとより少ない変数でエラーになる可能性があります．

おそらくGröbner基底簡約まわりで死んでるので，根本的に内部で同じことをやっているz3ならいけるかもしれません[^2].


# III. GF(p)の場合
簡単に言うとすべて$\pmod p$の世界です．ただし$p$は素数です．

基本的には II. と同じで，必要な部分に`GF(p)`を付け加えるだけでOK．

## III-a. 行列形式で表記できる場合

$$
\begin{cases}
x + 2y + 3z \equiv 0 \bmod 5 \\\\\
3x + 2y + z + 4 \equiv 0 \bmod 5 \\\\
x + y + 2z + 1 \equiv 0 \bmod 5
\end{cases}
$$


解は $(x, y, z) = (3, 1, 0)$です．

```python
A = Matrix(GF(5), [[1, 2, 3],
            [3, 2, 1],
            [1, 1, 2]])
v = vector(GF(5), [0, -4, -1])

print(A.solve_right(v))
```

## III-c. 行列形式では表記できない場合(変数がそこそこ多い場合)

III-b. は完全にこちらに吸収できるので省略．

また III-d. も同様のエラー持ちですが，z3で多項式環を処理する方法を筆者が知らない[^1]ため割愛．


$$
\begin{cases}
xy + z \equiv 0 \bmod 5 \\\\
x + yz \equiv 0 \bmod 5  \\\\
x + y + z + 1 \equiv 0 \bmod 5 
\end{cases}
$$

`QQ`を`GF(p)`に変えるだけでOK．

なお，$p=2$の場合にのみ，`PR = BooleanPolynomialRing(3, "x")`としたほうが早く動作します．

```python
PR = PolynomialRing(GF(5), 3, "x")
xvars = PR.gens()

eq1 = xvars[0] * xvars[1] + xvars[2]
eq2 = xvars[0] + xvars[1] * xvars[2]
eq3 = xvars[0] + xvars[1] + xvars[2] + 1

J = Ideal([eq1, eq2, eq3])
print(J.variety())
```






[^1]: 条件式に全て` % p`をつけて擬似的に$\bmod p$の世界を作り出し，かつ `0 <= x < p` などと条件追加すれば問題はないですが，それならばSageMathの多項式環クラス使ったほうが楽なのでは？と思っています．なので，少なくとも$p > 2$に関してはGF(p)とz3の親和性は高くないと考えています．

[^2]: z3はSageMathとは独立でC++実装だった気がするので同じエラーは吐かないと踏んでいます．
(つまりSageMath特有の`<class　'sage.rings.polynomial.multi_polynomial_sequence.PolynomialSequence_gf2'>`みたいな長文クラスは存在しない)
ただし，z3は構文解析などからやっている可能性があるのでSageMathより時間がかかるかもです．