-数学- はじめてSympyを使ってみた！　ーガロア群の計算を題材にしてー

f:id:jurupapa:20210506120723j:plain

この記事ではPython上で動作する数式処理システムSympyを使ってみた感想を書いていきます。

世の中ではPythonというプログラミング言語が大流行しています。スクリプト言語なので、簡単にプログラムを書いて試すことができること、データ処理、機械学習などのライブラリが充実していること、Jupyter notebookなど手軽に試せる環境が揃っていることなどが理由と思われます。

Pythonをベースにした数式処理システムはSympyとSageが有名です。SympyはMaximaなどと同じような比較的汎用な機能を提供しています。Pythonのライブラリの位置づけなので、他のPythonのライブラリと組み合わせて使うことも可能です。Sageは、専門的な機能(数論とか代数幾何とか)を提供する様々な数学ソフトをPythonから使えるようにしているため、ちゃんと使うには高度な数学を理解している必要があります。

普段はMaximaでの数式処理を紹介しているこのブログですが、今回は連休中に勉強したSympyについて使ってみた感想を書いていきます。題材としては「与えられた方程式のガロア群を具体的に計算する」プログラムを、Sympyで書き直してみました。

以前書いた上記の記事のシリーズで取り上げたアルゴリズムの流れです。ただ対称式の計算は時間がかかるので、ehitoさんに教えてもらった数値計算手法とグレブナー基底を使っています。

Maximaは数式処理言語ですが、Sympyは数式処理ライブラリです。この違いが最も大きく関係するのは数式の中の文字(x, y, a, b, cなど)の扱いです。Maximaでは言語の中の変数と数式の中の文字に区別はありません。一方Sympyではまず全く別物として扱い、出来るだけその区別をしなくても良い工夫をしています。

Sympyでは数式中で使う文字は全てオブジェクトです。そのため文字オブジェクトを生成して使います。

a=symbol('a')

右辺はシンボル型のオブジェクトで'a'という名前を持つ、数式の中で使える文字オブジェクトを生成しています。それを左辺の変数aに代入しておくことで、変数aを数式の中で使うと、文字オブジェクトsymbol('a')が実際には使われる、という仕組みです。ここは少し慣れが必要ですが、慣れればあまり苦にはなりませんでした。

Sympyにはかなり充実した置換群と置換が実装されています。今回の題材では、galResolvArray()という関数の定義で、n次対称群を生成し、その中の全ての置換で、解を表す変数リストを並び替える、という処理をほぼそのまま自然に描くことができています。symmetric(n)がn次対称群を生成する関数、gがその中の１つの要素で置換であり、g(varList)でvarListの要素をgで置換します。

for g in symmetric(len(varList) ):
VList.insert(0,phi(g(varList) ) )

また多項式を特別に定義して、多項式専用の操作をメソッドを使って行えます。係数のリストを取ったり、多項式の次数を求めることができます。この辺はよく整理されていて、使い慣れると便利です。

文字式の展開や因数分解などはほぼMaximaと同じやり方で実行できます。

Pythonを触ってみると、リスト処理の書き方が洗練されていることがわかります。この辺は慣れるとかなり便利です。例えば最後の方で定義しているmake_galois_group()という関数の中で、単拡大の原始元の最小多項式p4の解を作り出す置換を全て求めて１つのリストにまとめています。内包的と呼ばれる記法で書くと分かり易い記述が出来ます。

置換群が実装されていることのメリットですが、最小多項式を満たす可能な解の置換を全て集めて、１つの置換群を作り出すことが出来ます。このやり方で非常に自然にガロア群を置換群として生成し、求めることができます。

置換群には色々な性質が実装されており、例えば、is_solvableというメソッドを使うと、可解群かどうかを簡単に調べることができます。

多倍長数値計算についても多項式の全ての解を求めるnroots()メソッドや式への代入操作など、他のライブラリを導入することなくここで必要な計算を進めることができました。

用途を選ぶかもしれませんが、Sympyも結構使える印象を持ちました。

今回の記事を書くにあたって、日本語でのまとまったSympyの情報を以下のブログで勉強することができました。pianofisicaさんに感謝します。

pianofisica.hatenablog.com

以下は方程式のガロア群を計算するSympyプログラムです。f0にn次方程式を設定すれば、ggに置換群がもとまります。2〜4次方程式は早いのですが、5次方程式は簡単なもの(f0=x**5-2など)では動作しましたが、10分程度は時間がかかりました。

from sympy import *
from sympy.abc import a, b, c, d, e, f, x, y, z, V, k
from sympy.combinatorics import *
s1,s2,s3,s4,s5,s6=symbols('s1 s2 s3 s4 s5 s6')
SYMVARLIST=[s1,s2,s3,s4,s5,s6]
SOLVARLIST=[a,b,c,d,e,f]

f0=x**4 + x**3 + x**2 + x + 1
gbpoly_list=[ ]
factor(f0)

$$ x^{4} + x^{3} + x^{2} + x + 1$$

def symcoeff (f,var):
asollist=Array(SOLVARLIST[:poly(f).degree()])
p1=expand(Product(var-asollist[k],(k,0,len(asollist)-1) ).doit() )
return poly(f-p1,gens=[x]).all_coeffs()

def symvallist (f):
fpoly=poly(f)
symlist=SYMVARLIST[fpoly.degree()-1::-1]
res=[ ]
for deg in range(fpoly.degree() ):
res.append( (symlist[deg] , f.coeff(x,deg)*(-1)**(fpoly.degree()+deg) ) )
return res

symvallist(f0)

[(s4, 1), (s3, -1), (s2, 1), (s1, -1)]

gbpoly_list.extend(symcoeff(f0,x) )
Array(gbpoly_list)

$$ \left[a + b + c + d + 1,\, - a b - a c - a d - b c - b d - c d + 1, \, \\ a b c + a b d + a c d + b c d + 1,\, - a b c d + 1 \right]$$

def phi(args): return Sum(Array([1,-1,2,-5,0])[k]*Array(args)[k],(k,0,len(args)-1) ).doit()
def galResolvArray(varList):
VList=[ ]
GAList=[ ]
for g in symmetric(len(varList) ):
VList.insert(0,phi(g(varList) ) )
GAList.insert(0,g)
return( (Array(VList),Array(GAList) ) )

varray,galist=galResolvArray(SOLVARLIST[:poly(f0).degree()])
f2=Product(V-varray[k],(k,0,(varray.shape)[0]-1)).doit()

$$ \left(V - 2 a - b + c + 5 d\right) \left(V - 2 a - b + 5 c + d\right) \left(V - 2 a + b - c + 5 d\right) \left(V - 2 a + b + 5 c - d\right) \left(V - 2 a + 5 b - c + d\right) \left(V - 2 a + 5 b + c - d\right) \left(V - a - 2 b + c + 5 d\right) \left(V - a - 2 b + 5 c + d\right) \left(V - a + b - 2 c + 5 d\right) \left(V - a + b + 5 c - 2 d\right) \left(V - a + 5 b - 2 c + d\right) \left(V - a + 5 b + c - 2 d\right) \left(V + a - 2 b - c + 5 d\right) \left(V + a - 2 b + 5 c - d\right) \left(V + a - b - 2 c + 5 d\right) \left(V + a - b + 5 c - 2 d\right) \left(V + a + 5 b - 2 c - d\right) \left(V + a + 5 b - c - 2 d\right) \left(V + 5 a - 2 b - c + d\right) \left(V + 5 a - 2 b + c - d\right) \left(V + 5 a - b - 2 c + d\right) \left(V + 5 a - b + c - 2 d\right) \left(V + 5 a + b - 2 c - d\right) \left(V + 5 a + b - c - 2 d\right)$$

numsollist=list(zip(SOLVARLIST[:poly(f0).degree():],poly(f0).nroots(n=300)))
zcoeffs=poly(f2.subs(numsollist)).all_coeffs()
p4=factor(Poly(list(map(lambda x:x.round(),zcoeffs)),gens=V).expr).as_ordered_factors()[0]
gbpoly_list.append(p4)
p4

$$ V^{4} - 3 V^{3} - 26 V^{2} - 132 V + 1331$$

gbpoly_list.append(phi(SOLVARLIST[:poly(f0).degree()])-V)

Array(gbpoly_list)

$$ \left[a + b + c + d + 1,\, - a b - a c - a d - b c - b d - c d + 1 ,\, a b c + a b d + a c d + b c d + 1 ,\, - a b c d + 1 ,\, V^{4} - 3 V^{3} - 26 V^{2} - 132 V + 1331 ,\, - V + a - b + 2 c - 5 d \right]$$

gb_res=groebner( gbpoly_list, *SOLVARLIST[:poly(f0).degree()],V ,order='lex')
gb_res

$$ GroebnerBasis\left(\left( - 103 V^{3} - 1000 V^{2} - 512 V + 32681 a + 40722, \ - 199 V^{3} + 289 V^{2} + 2501 V + 32681 b + 34573, \ 202 V^{3} + 692 V^{2} - 5659 V + 32681 c - 31317, \ 100 V^{3} + 19 V^{2} + 3670 V + 32681 d - 11297, \ V^{4} - 3 V^{3} - 26 V^{2} - 132 V + 1331\right), \left( a, \ b, \ c, \ d, \ V\right)\right)$$

solVlist=solve(gb_res,SOLVARLIST[:poly(f0).degree()])

Dict(solVlist)

$$ \left\{ a : \frac{103 V^{3}}{32681} + \frac{1000 V^{2}}{32681} + \frac{512 V}{32681} - \frac{3702}{2971}, \ b : \frac{199 V^{3}}{32681} - \frac{289 V^{2}}{32681} - \frac{2501 V}{32681} - \frac{3143}{2971}, \ c : - \frac{202 V^{3}}{32681} - \frac{692 V^{2}}{32681} + \frac{5659 V}{32681} + \frac{2847}{2971}, \ d : - \frac{100 V^{3}}{32681} - \frac{19 V^{2}}{32681} - \frac{3670 V}{32681} + \frac{1027}{2971}\right\}$$

permV=list(map(lambda xx:xx.subs(solVlist),list(varray)))
Array(permV)

$$ \left[\begin{array}- \frac{15 V^{3}}{32681} - \frac{4905 V^{2}}{32681} - \frac{16891 V}{32681} + \frac{10404}{2971}, &\\ - \frac{687 V^{3}}{32681} + \frac{4118 V^{2}}{32681} + \frac{4200 V}{32681} + \frac{6491}{2971}, &\\ - \frac{1218 V^{3}}{32681} - \frac{6114 V^{2}}{32681} + \frac{7589 V}{32681} + \frac{28374}{2971}, &\\ \frac{917 V^{3}}{32681} + \frac{5730 V^{2}}{32681} - \frac{28440 V}{32681} - \frac{17469}{2971}, &\\ - \frac{1602 V^{3}}{32681} - \frac{958 V^{2}}{32681} + \frac{19641 V}{32681} + \frac{26138}{2971}, &\\ \frac{1205 V^{3}}{32681} + \frac{1863 V^{2}}{32681} - \frac{37479 V}{32681} - \frac{15792}{2971}, &\\ - \frac{219 V^{3}}{32681} - \frac{6251 V^{2}}{32681} + \frac{1767 V}{32681} + \frac{14044}{2971}, &\\ - \frac{81 V^{3}}{2971} + \frac{252 V^{2}}{2971} + \frac{2078 V}{2971} + \frac{10131}{2971}, &\\ - \frac{1116 V^{3}}{32681} - \frac{5441 V^{2}}{32681} - \frac{1740 V}{32681} + \frac{26554}{2971}, &\\ \frac{305 V^{3}}{32681} + \frac{1692 V^{2}}{32681} + \frac{27534 V}{32681} - \frac{6549}{2971}, &\\ - \frac{1500 V^{3}}{32681} - \frac{285 V^{2}}{32681} + \frac{10312 V}{32681} + \frac{24318}{2971}, &\\ \frac{593 V^{3}}{32681} - \frac{2175 V^{2}}{32681} + \frac{18495 V}{32681} - \frac{4872}{2971}, &\\ - \frac{620 V^{3}}{32681} - \frac{6654 V^{2}}{32681} + \frac{9927 V}{32681} + \frac{20034}{2971}, &\\ \frac{1515 V^{3}}{32681} + \frac{5190 V^{2}}{32681} - \frac{26102 V}{32681} - \frac{25809}{2971}, &\\ - \frac{314 V^{3}}{32681} - \frac{4635 V^{2}}{32681} - \frac{18060 V}{32681} + \frac{14574}{2971}, &\\ \frac{1107 V^{3}}{32681} + \frac{2498 V^{2}}{32681} + \frac{11214 V}{32681} - \frac{18529}{2971}, &\\ \frac{906 V^{3}}{32681} + \frac{2133 V^{2}}{32681} - \frac{38648 V}{32681} - \frac{11622}{2971}, &\\ \frac{192 V^{3}}{32681} - \frac{2578 V^{2}}{32681} + \frac{26655 V}{32681} + \frac{1118}{2971}, &\\ - \frac{1196 V^{3}}{32681} + \frac{1080 V^{2}}{32681} + \frac{28005 V}{32681} + \frac{16680}{2971}, &\\ \frac{1611 V^{3}}{32681} + \frac{3901 V^{2}}{32681} - \frac{29115 V}{32681} - \frac{25250}{2971} ,&\\ - \frac{890 V^{3}}{32681} + \frac{3099 V^{2}}{32681} + \frac{18 V}{32681} + \frac{11220}{2971}, &\\ \frac{1203 V^{3}}{32681} + \frac{1209 V^{2}}{32681} + \frac{8201 V}{32681} - \frac{17970}{2971}, &\\ \frac{714 V^{3}}{32681} + \frac{4711 V^{2}}{32681} - \frac{32622 V}{32681} - \frac{12740}{2971}, &\\ V\end{array}\right]$$

def make_part_list(lis,index_list):
return list(map(lambda index:lis[index],index_list) )

def make_galois_group(permV):
indexlist=[index for index,perm in enumerate(permV) if rem( expand(p4.subs(V,perm) ), p4)==0]
print(make_part_list(galist,indexlist) )
return PermutationGroup(list(map(Permutation,make_part_list(galist,indexlist) ) ) )

gg=make_galois_group(permV)
print(gg)