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

I. はじめに

corCTF 2022の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)$です．

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$ として扱われるので省略しても問題ないです．

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}
$$

を流用します．

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. 通りに実装すると以下のようになりますがエラーを吐きます．

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ならいけるかもしれません¹.

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)$です．

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で多項式環を処理する方法を筆者が知らない²ため割愛．

$$
\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")としたほうが早く動作します．

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

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