インターフェイスについて#

Sageの最大の特長は,多種多様なコンピュータ代数システムを,共通インターフェイスと本格的なプログラミング言語Pythonを用いて一つ屋根の下にまとめ上げている点だ. これによりSageでは由来の異なるオブジェクト群を組み合わせた演算処理が可能になっている.

ある1つのインターフェイスに対し, console メソッドと interact メソッドは,全く違った機能を果たしている. GAPを例にとって具体的に見てみよう:

  1. gap.console(): このメソッドはGAPのコンソールを起動し,命令実行の主体をGAP上へ移す. Sageの役割は,Linuxのbashシェルのような便利なプログラムランチャとしての役割に限定される.

  2. gap.interact(): このメソッドは,Sageオブジェクト群を使って動作しているGAPインスタンスと情報交換するための便利な径路を提供する. これを使うと,GAPセッション中に(対話型インタフェイスからでも)Sageオブジェクトをインポートすることができる.

GP/PARI#

数論関係の演算処理を主目的とするPARIは,コンパクトで非常に練れたプログラムで,高度に最適化されている. SageからPARIを使うには,大きく異なる2種類のインターフェイスを選ぶことができる:

  • gp -- "G o P ARI" インタープリタ

  • pari -- PARI Cライブラリ

以下の例では,同一の計算を二通りのやり方で実行している. 一見同じように見えても実は出力内容は同一ではないし,画面の奥で実際に行なわれている処理過程は二つの方法で全くと言っていいほど異なっているのだ.

sage: gp('znprimroot(10007)')
Mod(5, 10007)
sage: pari('znprimroot(10007)')
Mod(5, 10007)
>>> from sage.all import *
>>> gp('znprimroot(10007)')
Mod(5, 10007)
>>> pari('znprimroot(10007)')
Mod(5, 10007)

第一のやり方では,GPインタープリタが独立したサーバプロセスとして起動され,そのプロセスに文字列 'znprimroot(10007)' が送られる. GPは送られて来た文字列を評価し,結果を変数に格納する(変数は子GPプロセス配下の,開放されないメモリ空間内に確保される). その変数の値が画面表示されて仕上がりになる. 第二のやり方では、独立したプロセスが起動されることはなく,文字列 'znprimroot(10007)' はPARI Cライブラリ関数群によって処理される. 処理結果はPythonの確保しているヒープ上に配置され,その変数が参照されなくなると使っていたヒープメモリは開放される. 二つのやり方では、生成されるオブジェクトの型からして違っている:

sage: type(gp('znprimroot(10007)'))
<class 'sage.interfaces.gp.GpElement'>
sage: type(pari('znprimroot(10007)'))
<class 'cypari2.gen.Gen'>
>>> from sage.all import *
>>> type(gp('znprimroot(10007)'))
<class 'sage.interfaces.gp.GpElement'>
>>> type(pari('znprimroot(10007)'))
<class 'cypari2.gen.Gen'>

では,どちらの方法を選ぶべきだろうか? 答は目的による、としか言えない. GPインターフェイスはGP/PARIプログラムをそのまま動作させるのだから,GP/PARIのコマンド入力行から可能なことは全て出来る. 複雑なPARIプログラムを読み込んで走らせたい場合などには向いているだろう. これと比較すると,(Cライブラリを経由する)PARIインターフェイスはかなり制限がきつい. まだライブラリのメンバ関数の全てが実装されているわけではないし, 数値積分などを含むコードの多くがPARIインターフェイスからは使えない状態だ. 一方,PARIインターフェイスはGPインターフェイスより大幅に高速で頑強でもある.

(入力行を評価中にメモリ不足に陥った場合,GPインターフェイスは特に警告することなく自動的にスタックサイズを2倍して評価を再試行する. 必要とされるメモリ量を正しく見積っていなかったとしても処理が頓挫することはまずない. この有難い仕掛けは,標準のGPインタープリタには備えられていないようだ. PARI Cライブラリ インターフェイスについて言うと,こちらは生成したオジェクトを直ちにコピーしてPARIスタックから送り出してしまうので,スタックが積み上がることはない. しかし,どんなオブジェクトも大きさが100MBを越えてはならず,越えて生成された場合はスタックオーバーフローが起きる. このオブジェクトのコピーが,わずかながら実行効率上の不利を招いているのも確かである.)

まとめると,SageはPARI Cライブラリを利用してGP/PARIインタープリタと同様の機能を提供しているが,優秀なメモリ管理とプログラミング言語Pythonの援用という利点がある,ということになる.

ここで,PythonのリストからPARIのリストを作ってみよう.

sage: v = pari([1,2,3,4,5])
sage: v
[1, 2, 3, 4, 5]
sage: type(v)
<class 'cypari2.gen.Gen'>
>>> from sage.all import *
>>> v = pari([Integer(1),Integer(2),Integer(3),Integer(4),Integer(5)])
>>> v
[1, 2, 3, 4, 5]
>>> type(v)
<class 'cypari2.gen.Gen'>

PARIのオブジェクトは全て Gen 型になる. 各オブジェクトに埋め込まれているPARI由来のデータ型を取得するには,メンバ関数 type を使う.

sage: v.type()
't_VEC'
>>> from sage.all import *
>>> v.type()
't_VEC'

PARIで楕円曲線を生成するには, ellinit([1,2,3,4,5]) と入力する. Sageの方法も似ているが, ellinit を先の t_VEC \(v\) のような任意のPARIオブジェクトに対して呼び出すことができる点が違っている.

sage: e = v.ellinit()
sage: e.type()
't_VEC'
sage: pari(e)[:13]
[1, 2, 3, 4, 5, 9, 11, 29, 35, -183, -3429, -10351, 6128487/10351]
>>> from sage.all import *
>>> e = v.ellinit()
>>> e.type()
't_VEC'
>>> pari(e)[:Integer(13)]
[1, 2, 3, 4, 5, 9, 11, 29, 35, -183, -3429, -10351, 6128487/10351]

楕円曲線オブジェクトが生成できたので,これを使った計算を試みよう.

sage: e.elltors()
[1, [], []]
sage: e.ellglobalred()
[10351, [1, -1, 0, -1], 1, [11, 1; 941, 1], [[1, 5, 0, 1], [1, 5, 0, 1]]]
sage: f = e.ellchangecurve([1,-1,0,-1])
sage: f[:5]
[1, -1, 0, 4, 3]
>>> from sage.all import *
>>> e.elltors()
[1, [], []]
>>> e.ellglobalred()
[10351, [1, -1, 0, -1], 1, [11, 1; 941, 1], [[1, 5, 0, 1], [1, 5, 0, 1]]]
>>> f = e.ellchangecurve([Integer(1),-Integer(1),Integer(0),-Integer(1)])
>>> f[:Integer(5)]
[1, -1, 0, 4, 3]

GAP#

Sageには GAP 4.4.10が付属しており,群論を始めとする計算離散数学に対応している.

ここでは、例としてGAPの IdGroup 関数を取り上げることにしよう. この機能を使うには群論関係の小規模なデータベースが必要だが,標準ではインストールされない. これは別途インストールしておかなければならないから,後で手順を説明する.

sage: G = gap('Group((1,2,3)(4,5), (3,4))')
sage: G
Group( [ (1,2,3)(4,5), (3,4) ] )
sage: G.Center()
Group( () )
sage: G.IdGroup()
[ 120, 34 ]
sage: G.Order()
120
>>> from sage.all import *
>>> G = gap('Group((1,2,3)(4,5), (3,4))')
>>> G
Group( [ (1,2,3)(4,5), (3,4) ] )
>>> G.Center()
Group( () )
>>> G.IdGroup()
[ 120, 34 ]
>>> G.Order()
120

上と同じ処理を,GAPインタ=フェイスを明示的には呼び出さずにSageから実行するには:

sage: G = PermutationGroup([[(1,2,3),(4,5)],[(3,4)]])
sage: G.center()
Subgroup generated by [()] of (Permutation Group with generators [(3,4), (1,2,3)(4,5)])
sage: G.group_id()
[120, 34]
sage: n = G.order(); n
120
>>> from sage.all import *
>>> G = PermutationGroup([[(Integer(1),Integer(2),Integer(3)),(Integer(4),Integer(5))],[(Integer(3),Integer(4))]])
>>> G.center()
Subgroup generated by [()] of (Permutation Group with generators [(3,4), (1,2,3)(4,5)])
>>> G.group_id()
[120, 34]
>>> n = G.order(); n
120

(GAPの機能の一部は,二種類の標準外Sageパッケージをインストールしないと使うことができない. コマンド行から sage -optional と入力すると,インストール可能なパッケージの一覧が表示される. その一覧に gap\_packages-x.y.z といった項目があるから, sage -i gap\_packages-x.y.z を実行してパッケージをインストールする. 一部のGPLではないGAPパッケージは,GAPウェブサイト [GAPkg] からダウンロードし, $SAGE_ROOT/local/lib/gap-4.4.10/pkg で解凍してからインストールする必要がある.)

Singular#

Singularは,グレブナー基底,多変数多項式のgcd,平面曲線のRieman-Roch空間に対する基底、因数分解などを始めとする各種処理のための,大規模で十分に枯れたライブラリを提供する. 実例として,多変数多項式の因数分解をSageからSingularへのインターフェイスを使って実行してみよう(.... は入力しないこと):

sage: R1 = singular.ring(0, '(x,y)', 'dp')
sage: R1
polynomial ring, over a field, global ordering
//   coefficients: QQ
//   number of vars : 2
//        block   1 : ordering dp
//                  : names    x y
//        block   2 : ordering C
sage: f = singular('9*y^8 - 9*x^2*y^7 - 18*x^3*y^6 - 18*x^5*y^6 +'
....: '9*x^6*y^4 + 18*x^7*y^5 + 36*x^8*y^4 + 9*x^10*y^4 - 18*x^11*y^2 -'
....: '9*x^12*y^3 - 18*x^13*y^2 + 9*x^16')
>>> from sage.all import *
>>> R1 = singular.ring(Integer(0), '(x,y)', 'dp')
>>> R1
polynomial ring, over a field, global ordering
//   coefficients: QQ
//   number of vars : 2
//        block   1 : ordering dp
//                  : names    x y
//        block   2 : ordering C
>>> f = singular('9*y^8 - 9*x^2*y^7 - 18*x^3*y^6 - 18*x^5*y^6 +'
... '9*x^6*y^4 + 18*x^7*y^5 + 36*x^8*y^4 + 9*x^10*y^4 - 18*x^11*y^2 -'
... '9*x^12*y^3 - 18*x^13*y^2 + 9*x^16')

\(f\) を定義できたので,その内容を表示してから因数分解を試みる.

sage: f
9*x^16-18*x^13*y^2-9*x^12*y^3+9*x^10*y^4-18*x^11*y^2+36*x^8*y^4+18*x^7*y^5-18*x^5*y^6+9*x^6*y^4-18*x^3*y^6-9*x^2*y^7+9*y^8
sage: f.parent()
Singular
sage: F = f.factorize(); F
[1]:
   _[1]=9
   _[2]=x^6-2*x^3*y^2-x^2*y^3+y^4
   _[3]=-x^5+y^2
[2]:
   1,1,2
sage: F[1][2]
x^6-2*x^3*y^2-x^2*y^3+y^4
>>> from sage.all import *
>>> f
9*x^16-18*x^13*y^2-9*x^12*y^3+9*x^10*y^4-18*x^11*y^2+36*x^8*y^4+18*x^7*y^5-18*x^5*y^6+9*x^6*y^4-18*x^3*y^6-9*x^2*y^7+9*y^8
>>> f.parent()
Singular
>>> F = f.factorize(); F
[1]:
   _[1]=9
   _[2]=x^6-2*x^3*y^2-x^2*y^3+y^4
   _[3]=-x^5+y^2
[2]:
   1,1,2
>>> F[Integer(1)][Integer(2)]
x^6-2*x^3*y^2-x^2*y^3+y^4

GAP 節におけるGAPの実行例のように,Singularインターフェイスを直には使わず上の因数分解を行なうこともできる(Sageが実際の計算に裏でSingularインターフェイスを使っていることに変わりない). 以下の例でも, .... は入力しないこと:

sage: x, y = QQ['x, y'].gens()
sage: f = (9*y^8 - 9*x^2*y^7 - 18*x^3*y^6 - 18*x^5*y^6 + 9*x^6*y^4
....: + 18*x^7*y^5 + 36*x^8*y^4 + 9*x^10*y^4 - 18*x^11*y^2 - 9*x^12*y^3
....: - 18*x^13*y^2 + 9*x^16)
sage: factor(f)
(9) * (-x^5 + y^2)^2 * (x^6 - 2*x^3*y^2 - x^2*y^3 + y^4)
>>> from sage.all import *
>>> x, y = QQ['x, y'].gens()
>>> f = (Integer(9)*y**Integer(8) - Integer(9)*x**Integer(2)*y**Integer(7) - Integer(18)*x**Integer(3)*y**Integer(6) - Integer(18)*x**Integer(5)*y**Integer(6) + Integer(9)*x**Integer(6)*y**Integer(4)
... + Integer(18)*x**Integer(7)*y**Integer(5) + Integer(36)*x**Integer(8)*y**Integer(4) + Integer(9)*x**Integer(10)*y**Integer(4) - Integer(18)*x**Integer(11)*y**Integer(2) - Integer(9)*x**Integer(12)*y**Integer(3)
... - Integer(18)*x**Integer(13)*y**Integer(2) + Integer(9)*x**Integer(16))
>>> factor(f)
(9) * (-x^5 + y^2)^2 * (x^6 - 2*x^3*y^2 - x^2*y^3 + y^4)

Maxima#

Maximaは,LISP言語の実装の一種と共にSageに同梱されてくる. (Maximaが標準でプロットに用いる)gnuplotパッケージは,Sageでも非標準パッケージとして公開されている. Maximaが得意とするのは,記号処理である.Maximaは関数の微分と積分を記号的に実行し,1階の常微分方程式(ODE)と大半の線形2次ODEを解くことができるし,任意次数の線形ODEをラプラス変換することもできる. さらにMaximaには多様な特殊関数も組込まれており,gnuplotを介したプロット機能も備えている. (掃き出し法や固有値問題などに始まる)行列操作や行列方程式の解法を実行し,多項式方程式を解くことも可能だ.

Sage/Maximaインターフェイスの使い方を例示するため,ここでは \(i,j\) 要素の値が \(i/j\) で与えられる行列を作成してみよう. ただし \(i,j=1,\ldots,4\) とする.

sage: f = maxima.eval('ij_entry[i,j] := i/j')
sage: A = maxima('genmatrix(ij_entry,4,4)'); A
matrix([1,1/2,1/3,1/4],[2,1,2/3,1/2],[3,3/2,1,3/4],[4,2,4/3,1])
sage: A.determinant()
0
sage: A.echelon()
matrix([1,1/2,1/3,1/4],[0,0,0,0],[0,0,0,0],[0,0,0,0])
sage: A.eigenvalues()
[[0,4],[3,1]]
sage: A.eigenvectors().sage()
[[[0, 4], [3, 1]], [[[1, 0, 0, -4], [0, 1, 0, -2], [0, 0, 1, -4/3]], [[1, 2, 3, 4]]]]
>>> from sage.all import *
>>> f = maxima.eval('ij_entry[i,j] := i/j')
>>> A = maxima('genmatrix(ij_entry,4,4)'); A
matrix([1,1/2,1/3,1/4],[2,1,2/3,1/2],[3,3/2,1,3/4],[4,2,4/3,1])
>>> A.determinant()
0
>>> A.echelon()
matrix([1,1/2,1/3,1/4],[0,0,0,0],[0,0,0,0],[0,0,0,0])
>>> A.eigenvalues()
[[0,4],[3,1]]
>>> A.eigenvectors().sage()
[[[0, 4], [3, 1]], [[[1, 0, 0, -4], [0, 1, 0, -2], [0, 0, 1, -4/3]], [[1, 2, 3, 4]]]]

使用例をもう一つ示す:

sage: A = maxima("matrix ([1, 0, 0], [1, -1, 0], [1, 3, -2])")
sage: eigA = A.eigenvectors()
sage: V = VectorSpace(QQ,3)
sage: eigA
[[[-2,-1,1],[1,1,1]],[[[0,0,1]],[[0,1,3]],[[1,1/2,5/6]]]]
sage: v1 = V(sage_eval(repr(eigA[1][0][0]))); lambda1 = eigA[0][0][0]
sage: v2 = V(sage_eval(repr(eigA[1][1][0]))); lambda2 = eigA[0][0][1]
sage: v3 = V(sage_eval(repr(eigA[1][2][0]))); lambda3 = eigA[0][0][2]

sage: M = MatrixSpace(QQ,3,3)
sage: AA = M([[1,0,0],[1, - 1,0],[1,3, - 2]])
sage: b1 = v1.base_ring()
sage: AA*v1 == b1(lambda1)*v1
True
sage: b2 = v2.base_ring()
sage: AA*v2 == b2(lambda2)*v2
True
sage: b3 = v3.base_ring()
sage: AA*v3 == b3(lambda3)*v3
True
>>> from sage.all import *
>>> A = maxima("matrix ([1, 0, 0], [1, -1, 0], [1, 3, -2])")
>>> eigA = A.eigenvectors()
>>> V = VectorSpace(QQ,Integer(3))
>>> eigA
[[[-2,-1,1],[1,1,1]],[[[0,0,1]],[[0,1,3]],[[1,1/2,5/6]]]]
>>> v1 = V(sage_eval(repr(eigA[Integer(1)][Integer(0)][Integer(0)]))); lambda1 = eigA[Integer(0)][Integer(0)][Integer(0)]
>>> v2 = V(sage_eval(repr(eigA[Integer(1)][Integer(1)][Integer(0)]))); lambda2 = eigA[Integer(0)][Integer(0)][Integer(1)]
>>> v3 = V(sage_eval(repr(eigA[Integer(1)][Integer(2)][Integer(0)]))); lambda3 = eigA[Integer(0)][Integer(0)][Integer(2)]

>>> M = MatrixSpace(QQ,Integer(3),Integer(3))
>>> AA = M([[Integer(1),Integer(0),Integer(0)],[Integer(1), - Integer(1),Integer(0)],[Integer(1),Integer(3), - Integer(2)]])
>>> b1 = v1.base_ring()
>>> AA*v1 == b1(lambda1)*v1
True
>>> b2 = v2.base_ring()
>>> AA*v2 == b2(lambda2)*v2
True
>>> b3 = v3.base_ring()
>>> AA*v3 == b3(lambda3)*v3
True

最後に,Sage経由で openmath を使ってプロットを行なう際の手順を紹介する. 以下の例題の多くは,Maximaのレファレンスマニュアルのものを修正したものだ.

関数の2次元プロットをするには( ... は入力しない)

sage: maxima.plot2d('[cos(7*x),cos(23*x)^4,sin(13*x)^3]','[x,0,1]',  # not tested
....:     '[plot_format,openmath]')
>>> from sage.all import *
>>> maxima.plot2d('[cos(7*x),cos(23*x)^4,sin(13*x)^3]','[x,0,1]',  # not tested
...     '[plot_format,openmath]')

次の「ライブ」3次元プロットは,マウスで動かすことができる( .... は入力しない):

sage: maxima.plot3d ("2^(-u^2 + v^2)", "[u, -3, 3]", "[v, -2, 2]",  # not tested
....:     '[plot_format, openmath]')
sage: maxima.plot3d("atan(-x^2 + y^3/4)", "[x, -4, 4]", "[y, -4, 4]",  # not tested
....:     "[grid, 50, 50]",'[plot_format, openmath]')
>>> from sage.all import *
>>> maxima.plot3d ("2^(-u^2 + v^2)", "[u, -3, 3]", "[v, -2, 2]",  # not tested
...     '[plot_format, openmath]')
>>> maxima.plot3d("atan(-x^2 + y^3/4)", "[x, -4, 4]", "[y, -4, 4]",  # not tested
...     "[grid, 50, 50]",'[plot_format, openmath]')

次に有名なメビウスの帯を3次元プロットしてみよう( .... は入力しない).

sage: maxima.plot3d("[cos(x)*(3 + y*cos(x/2)), sin(x)*(3 + y*cos(x/2)), y*sin(x/2)]",  # not tested
....:     "[x, -4, 4]", "[y, -4, 4]", '[plot_format, openmath]')
>>> from sage.all import *
>>> maxima.plot3d("[cos(x)*(3 + y*cos(x/2)), sin(x)*(3 + y*cos(x/2)), y*sin(x/2)]",  # not tested
...     "[x, -4, 4]", "[y, -4, 4]", '[plot_format, openmath]')

プロットの最後の例は,あの「クラインの壺」である( .... は入力しない):

sage: _ = maxima("expr_1: 5*cos(x)*(cos(x/2)*cos(y) + sin(x/2)*sin(2*y)+ 3.0) - 10.0")
sage: _ = maxima("expr_2: -5*sin(x)*(cos(x/2)*cos(y) + sin(x/2)*sin(2*y)+ 3.0)")
sage: _ = maxima("expr_3: 5*(-sin(x/2)*cos(y) + cos(x/2)*sin(2*y))")
sage: maxima.plot3d ("[expr_1, expr_2, expr_3]", "[x, -%pi, %pi]",  # not tested
....:     "[y, -%pi, %pi]", "['grid, 40, 40]", '[plot_format, openmath]')
>>> from sage.all import *
>>> _ = maxima("expr_1: 5*cos(x)*(cos(x/2)*cos(y) + sin(x/2)*sin(2*y)+ 3.0) - 10.0")
>>> _ = maxima("expr_2: -5*sin(x)*(cos(x/2)*cos(y) + sin(x/2)*sin(2*y)+ 3.0)")
>>> _ = maxima("expr_3: 5*(-sin(x/2)*cos(y) + cos(x/2)*sin(2*y))")
>>> maxima.plot3d ("[expr_1, expr_2, expr_3]", "[x, -%pi, %pi]",  # not tested
...     "[y, -%pi, %pi]", "['grid, 40, 40]", '[plot_format, openmath]')