Sharedresdiary / 20171211 / advent.ipynbOpen in CoCalc

Sagemathで数学

この記事は

この記事はアドベントカレンダー「数学とコンピュータ」の12月13日の記事です. 前の記事は @n_kats_ さん https://qiita.com/n_kats_/items/99eb63e2ede23e05e059 ,次の記事は @mod_poppo さんhttps://qiita.com/mod_poppo です.

フリーの総合数学ソフトSagemathの紹介と,ちょっとしたデモをご覧いただきます.

Sagemathとは

Sagemath http://www.sagemath.org/ は,「Mathematica, Maple, Magma, Matlabなどを置き換えうるフリー の総合数学ソフトウェア」を目指して開発が進められているフリーソフトです.

このSagemath開発プロジェクトはW. Steinさんによってはじめられ,現在も活発な開発が続いています.

Sagemathの実態は,Python2に様々な既存の数学ソフトウェア・もしくはSagemathのために新規に開発された数学ソフトウェアを統合したものです.C/Fortranのライブラリとして取り込まれたもの(例えばLapack/linalg/pari...)が多いですが,例えばpexpectを使って,仮想端末越しに操作している例もあります.

Sagemathを使ってみたい

Sagemathのバイナリは,各種Linuxディストリビューション,Macのほか,Windows版も提供されています.かなり巨大なのでダウンロードには少し待ち時間がかかるでしょう.(ダウンロードサイト:http://www.sagemath.org/download.html

自分のパソコンにインストールしなくても,Cocalcというサービス http://cocal.com を無料で使い始めることができます.このサービスを運営している sagemath.inc の代表は,Sagemathの開発者 W. Stein さんその人です.

Sagemathでなにができる?

Sagemathは,上述のように,著名な数学ソフトウェアを置き換えることを目標 に開発されています.ですから,そういったソフトウェアでできることはたい ていできます.

  • 数値計算
  • 数式処理(多項式や特殊関数の演算,式変形)
  • グラフィックス(グラフ,3d表示,レイトレーシング)

それだけではなく,Sagemathの開発者は数学,特に数論や組み合わせ論,幾何 学の研究者が多数参画しています.W. Steinさんが,モジュラー形式やアーベ ル多様体と呼ばれる対象の,特に数論的な性質を調べる分野(いわゆる数論幾 何)の専門家であることから,そういった方面に特化した計算が得意です. 例えば,

の計算が,「コマンド一発」で実行できます.

少し試してみましょう.

有限体

有限体の計算を,ElGamal暗号の素朴な実装をしながら見てみましょう.まず,pp2q+12q+1, qqも素数,の形に取りたいので,探してみます.

def search_sg_prime(f, t):
    """Search Sophie-German prime p of the form 2q+1, q is a prime."""
    for q in primes(f, t):
        p = 2*q + 1
        if is_prime(p):
            return(p, q)
    return(None)
p, q = search_sg_prime(10^4, 10^4+100); p,q
(20123, 10061)
def system_setup(p, q):
    F = FiniteField(p)
    gamma = F.multiplicative_generator() # 法pの原始根
    g = gamma^2 # p=2q+1から,2 = (p-1)/q, つまりgは位数q
    x = ZZ.random_element(0, q-1)
    h = g^x
    return(g, h, x)
g, h, x = system_setup(p, q); g, h, x
(4, 13145, 410)

すると,(G=g,q=10061,g,h)(G = \langle g\rangle, q=10061, g, h )が公開鍵,xxが秘密鍵.

暗号化:平文 mGm\in Gに対して r{0,,q1}r\in \{0,\dots, q-1\} をランダムに取って,(c1,c2)=(gr,mhr)(c_1, c_2) = (g^r, m h^r) で定まる (c1,c2)(c_1, c_2) が暗号文.

復号化:暗号文 (c1,c2)(c_1, c_2) に対して,m=c2(c1x)1m = c_2 (c_1^x)^{-1}.

m = g^1213; m # 12月13日なので.
16301
def our_encrypt(m, g, h, x):
    r = ZZ.random_element(0, q-1)
    c1 = g^r; c2 = m*h^r
    return(c1, c2)
c1, c2 = our_encrypt(m, g, h, x); c1, c2
(3016, 4406)
def our_decrypt(c1, c2, g, h, x):
    return(c2*c1^(-x))
our_decrypt(c1, c2, g, h, x)
16301

無事平文が得られていることがわかります.秘密鍵 xx が離散対数の計算で分かってしまうと台無しですが,このくらいの大きさでは簡単に計算できてしまいます:

h.log(g)
410

代数体

虚2次体の類数

虚2次体の類数の計算からはじめます.類数が9で割り切れるものを列挙してみましょう.そして,類数が9で割り切れるということは,イデアル類群が位数9の巡回群もしくは(3,3)型の有限アーベル群を部分群として含むと言うことですが,どちらなのかを見てみましょう.

実は,虚2次体のイデアル類群のpp Sylow群は,同じ位数なら巡回群の方がより頻繁に起きるという経験則がしられています.これを定式化したものがCohen-Lenstra予想で,現在でも未解決です.

判別式の絶対値が500以下の虚2次体で,類数が9で割り切れるものを列挙.

for d in range(1, 500):
    if (is_fundamental_discriminant(-d)):
        h = QuadraticField(-d).class_number()
        if (h%9 == 0):
            print(-d, h)
(-199, 9) (-335, 18) (-367, 9) (-419, 9) (-491, 9)

上で計算した類数が9で割り切れる虚2次体のイデアル類群を列挙すると,どれも位数9の巡回群.

[(-d, QuadraticField(-d).class_group()) for d in [199, 335, 367, 419, 491]]
[(-199, Class group of order 9 with structure C9 of Number Field in a with defining polynomial x^2 + 199), (-335, Class group of order 18 with structure C18 of Number Field in a with defining polynomial x^2 + 335), (-367, Class group of order 9 with structure C9 of Number Field in a with defining polynomial x^2 + 367), (-419, Class group of order 9 with structure C9 of Number Field in a with defining polynomial x^2 + 419), (-491, Class group of order 9 with structure C9 of Number Field in a with defining polynomial x^2 + 491)]

さて,ちょっと目先を変えて,虚2次体の類数をグラフに描くと,何が見えてくるでしょうか?.

fund_ds = [d for d in range(1, 2000) if is_fundamental_discriminant(-d)];
pts = sum([point([d, QuadraticField(-d).class_number()]) for d in fund_ds])
pts.show()

円分体の類数の計算.

小さい素数分体で試すと問題ない:

[(p, CyclotomicField(p).class_number()) for p in primes(20)]
[(2, 1), (3, 1), (5, 1), (7, 1), (11, 1), (13, 1), (17, 1), (19, 1)]

少し pp が大きくなると,てきめんに結果が返ってこない.特に工夫せずに類数やイデアル類群が計算できるのは,有理数体上十数次の拡大体までです(経験的に):

# [(p, CyclotomicField(p).class_number()) for p in primes(20, 30)] # 結果は戻ってこない

円分体(いま考えるのは素数 p>2p>2 に対する素数分体)の類数 h(p)h(p) は,h(p)=h(p)+h(p)h(p) = h(p)^+ h(p)^- のように正整数の積であらわされることが知られています.h(p)+h(p)^+pp 分体の最大実部分体 Q(ζp+ζp1)\mathbb{Q}(\zeta_p+\zeta_p^{-1}) の類数.h(p)h(p)^- は上の式で定義され,pp分体の相対類数とも呼ばれます. 実は,円分体の相対類数には,次のような解析的な表示式が知られています: h(p)=Qpwpχ(12B1,χ), h(p)^- = Q_p w_p \prod_{\chi} \left(-\frac{1}{2} B_{1, \chi}\right), QpQ_pは単数指数(unit index)と呼ばれる量で,今の場合22, wpw_p は1の冪根の個数で,いまの場合はpp です. 積は法ppの奇なDirichlet指標に渡り,B1,χB_{1,\chi}χ\chi に関する1番目の一般化Bernoulli数.(この表式については,例えばL. C. Washigton, Introduction to Cyclotomic Fields, SpringerのTheorem 4.17をご覧下さい).

一方,最大実部分体の類数は大変難しく,最近までの結果をまとめると,,p151p \le 151 のときh(p)+=1h(p)^+=1. GRHを仮定すると,p263p\le 263, p163,191,229,257p\neq 163, 191, 229, 257のときhp+=1h_p^+=1, h163+=4h_{163}^+=4, h191+=11h_{191}^+=11, h229+=3h_{229}^+=3, h257+=3h_{257}^+=3 が知られています(cf. J. C. Miller, RIMS Kokyuroku Bessatsu B64, 2017).

これらをまとめて,次のような素数分体の類数を計算する関数が書けるでしょう:

def hpminus(p):
    rv = prod(-(1/2)*chi.bernoulli(1) for chi in DirichletGroup(p) if chi.is_odd())
    return(ZZ(rv*p*2))
def hpplus(p):
    """Without GRH."""
    if p <= 151:
        return(1)
    return(None)
def hp(p):
    """class number of the p-th cyclotomic field."""
    return(hpplus(p) * hpminus(p))

さきほど結果がかえってこないと言った類数の計算も,解析的類数公式を使えば可能です.相対類数は素早く増加することが見て取れます.

[(p, hp(p)) for p in primes(20, 100)]
[(23, 3), (29, 8), (31, 9), (37, 37), (41, 121), (43, 211), (47, 695), (53, 4889), (59, 41241), (61, 76301), (67, 853513), (71, 3882809), (73, 11957417), (79, 100146415), (83, 838216959), (89, 13379363737), (97, 411322824001)]

h(p)h(p)^-の表は,Washingtonの上掲書Tables §3にもあります.

有限体上の楕円曲線

有限体上の楕円曲線のFpF_p有理点群が巡回群になるものを探し,有限体の乗法群とときと同じように,Elgamal暗号を構成してみましょう. FpF_p 上の楕円曲線 EEFpF_p有理点群 E(Fp)E(F_p) は巡回群であるか,正整数m1,m2m_1, m_2m1m2m_1 \mid m_2, m1gcd(#E(Fp),p1)m_1 \mid \gcd(\# E(F_p), p-1) となることが知られています(Cassels).したがって特に,位数が平方因子を含まなければ巡回群です.そのようなものを探索します.

p = next_prime(10^4); Fp = FiniteField(p); Fp
Finite Field of size 10007
def find_ec(p):
    F = FiniteField(p)
    while True:
        a = F(ZZ.random_element(0, p))
        b = F(ZZ.random_element(0, p))
        try:
            E = EllipticCurve([a, b])
            if E.order().is_squarefree():
                return(E)
            else:
                next
        except ArrithmeticError:
            next

あとは有限体のときと同様,システムのセットアップ,暗号化,復号化の手続きを定義します.

def system_setup_ell(p):
    E = find_ec(p)
    g = E.gen(0)
    x = ZZ.random_element(0, E.order()-1)
    h = x*g
    return(g, h, x, E)
def our_encrypt_ell(m, r, g, h, E):
    r = ZZ.random_element(0, E.order()-1)
    c1 = r*g; c2 = m + r*h
    return(c1, c2)
def our_decrypt_ell(c1, c2, x):
    return(c2 - x*c1)
g, h, x, E = system_setup_ell(next_prime(10^4))
m = 1213*g; m # やはり12月13日にちなんで.
(7357 : 8126 : 1)
r = ZZ.random_element(0, E.order()-1); r
5808
c1,c2 = our_encrypt_ell(m, r, g, h, E)
our_decrypt_ell(c1, c2, x)
(7357 : 8126 : 1)

有理数体上の楕円曲線

1213にちなんだ合同数問題を考えてみましょう.1213は合同数でしょうか?(正整数nnが合同数であるとは,直角三角形で,3辺の長さが有理数かつ,面積がnnであるようなものが存在することを言います.例えばWikipedia https://en.wikipedia.org/wiki/Congruent_number 参照).簡単な考察から,楕円曲線 E ⁣:y2=x3n2xE\colon y^2 = x^3 - n^2 xy0y\neq 0 である有理点を持てばよいと言うことが分かります.特に,EEの有理点群(Mordell-Weil群,これは有限生成アーベル群であることが知られている)が正のランクを持てば nn は合同数です.

n=1213n=1213 の場合を考えます:

E = EllipticCurve([-1213^2, 0]); E
Elliptic Curve defined by y^2 = x^3 - 1471369*x over Rational Field
E.torsion_points()
[(-1213 : 0 : 1), (0 : 0 : 1), (0 : 1 : 0), (1213 : 0 : 1)]

トーション点はすぐに求められましたが,有理点群(Mordell-Weil群)の生成元の計算には失敗してしまいます:

# E.gens() # 失敗する.メッセージに従って,2-descentのパラメータを調整して再挑戦
# E.gens(descent_second_limit=14) # 残念ながらこちらも失敗する.

analytic rankは計算でき,1になる.BSDが正しければ,EのMordell-Weil群はランク1で,従って n=1213n=1213 は合同数です.

E.analytic_rank()
1
E.analytic_rank_upper_bound()
1

Cremonaのmwrankを直接起動してみます.

print(E.mwrank())
Curve [0,0,0,-1471369,0] : 3 points of order 2: [-1213:0:1], [0:0:1], [1213:0:1] **************************** * Using 2-isogeny number 1 * **************************** Using 2-isogenous curve [0,7278,0,1471369,0] (minimal model [0,0,0,-16185059,24986788358]) ------------------------------------------------------- First step, determining 1st descent Selmer groups ------------------------------------------------------- After first local descent, rank bound = 1 rk(S^{phi}(E'))= 2 rk(S^{phi'}(E))= 1 ------------------------------------------------------- Second step, determining 2nd descent Selmer groups ------------------------------------------------------- After second local descent, rank bound = 1 rk(phi'(S^{2}(E)))= 2 rk(phi(S^{2}(E')))= 1 rk(S^{2}(E))= 3 rk(S^{2}(E'))= 2 **************************** * Using 2-isogeny number 2 * **************************** Using 2-isogenous curve [0,0,0,5885476,0] ------------------------------------------------------- First step, determining 1st descent Selmer groups ------------------------------------------------------- After first local descent, rank bound = 1 rk(S^{phi}(E'))= 2 rk(S^{phi'}(E))= 1 ------------------------------------------------------- Second step, determining 2nd descent Selmer groups ------------------------------------------------------- After second local descent, rank bound = 1 rk(phi'(S^{2}(E)))= 2 rk(phi(S^{2}(E')))= 1 rk(S^{2}(E))= 3 rk(S^{2}(E'))= 2 **************************** * Using 2-isogeny number 3 * **************************** Using 2-isogenous curve [0,-7278,0,1471369,0] (minimal model [0,0,0,-16185059,-24986788358]) ------------------------------------------------------- First step, determining 1st descent Selmer groups ------------------------------------------------------- After first local descent, rank bound = 1 rk(S^{phi}(E'))= 3 rk(S^{phi'}(E))= 0 ------------------------------------------------------- Second step, determining 2nd descent Selmer groups ------------------------------------------------------- After second local descent, rank bound = 1 rk(phi'(S^{2}(E)))= 3 rk(phi(S^{2}(E')))= 0 rk(S^{2}(E))= 3 rk(S^{2}(E'))= 2 After second local descent, combined upper bound on rank = 1 Third step, determining E(Q)/phi(E'(Q)) and E'(Q)/phi'(E(Q)) ------------------------------------------------------- 1. E(Q)/phi(E'(Q)) ------------------------------------------------------- (c,d) =(3639,2942738) (c',d')=(-7278,1471369) First stage (no second descent yet)... (-1,0,3639,0,-2942738): no rational point found (hlim=8) (2426,0,3639,0,1213): no rational point found (hlim=8) (-2,0,3639,0,-1471369): no rational point found (hlim=8) (1213,0,3639,0,2426): no rational point found (hlim=8) After first global descent, this component of the rank has lower bound 0 and upper bound 1 (difference = 1) Second descent will attempt to reduce this Second stage (using second descent)... d1=-1: Second descent inconclusive for d1=-1: ELS descendents exist but no rational point found d1=2426: Second descent inconclusive for d1=2426: ELS descendents exist but no rational point found d1=-2: Second descent inconclusive for d1=-2: ELS descendents exist but no rational point found d1=1213: Second descent inconclusive for d1=1213: ELS descendents exist but no rational point found After second global descent, this component of the rank has lower bound 0 and upper bound 1 (difference = 1) ------------------------------------------------------- 2. E'(Q)/phi'(E(Q)) ------------------------------------------------------- This component of the rank is 0 ------------------------------------------------------- Summary of results: ------------------------------------------------------- 0 <= rank(E) <= 1 #E(Q)/2E(Q) >= 4 Information on III(E/Q): #III(E/Q)[phi'] = 1 #III(E/Q)[2] <= 2 Information on III(E'/Q): #phi'(III(E/Q)[2]) <= 2 #III(E'/Q)[phi] <= 2 #III(E'/Q)[2] <= 2 Used descent via 2-isogeny with isogenous curve E' = [0,0,0,-16185059,-24986788358] 0 <= rank <= 1 Searching for points (bound = 8)...done: found points which generate a subgroup of rank 0 and regulator 1 Processing points found during 2-descent...done: now regulator = 1 Regulator = 1 The rank has not been completely determined, only a lower bound of 0 and an upper bound of 1. (0.48 seconds)
# print(E.mwrank('-b 15')) # パラメータを調整してCremonaのmwrankを起動しても,やはりうまくいかない.

実は,上のWikipediaの記事(に引用されている論文 https://eudml.org/doc/183796 )によれば,法8で5と合同な素数(1213はそうです)は合同数だそうです.

長くなったので,モジュラー形式については稿をあらためて書きたいと思います.(待ちきれない方は,thematic tutorialからどうぞ).


フリーソフトとしての数学ソフトウェア開発

フリーソフトウェアとしての数学ソフトウェアの開発にはいくつか困難があるように思います.

数学の非自明な計算をコンピュータに実装できる人材を安定して雇用するには,産学官の補助金などに頼るプロジェクトは脆弱すぎるのです.また,アマチュアが自分の興味に基づいて,断続的に貢献するというのでは,商用の数学ソフトウェアを置き換えるものには達しえないのではないか,という懸念がどうしてもおきるのです.そういった懸念を率直に語ったblog記事が W. Stein さんの http://sagemath.blogspot.jp/2014/08/you-dont-really-think-that-sage-has.html です.

Cocalcは,Sagemathを気軽に試せる環境を提供し,無料から使い始めて,必要に応じて少額のsubscription planを購入してもらい,利益を開発者の雇用に回す,という方針のようです. 個人的に応援しています.今回のadvent calendarの記事で,Sagemath, そしてCocalcに興味を持って頂けたら,ぜひ一度お試し下さい(もちろん私はなんの利害もありません).

文献など

必要に応じて引用しましたが,Sagemathについて日本語で読めるものとして

その他気づいたものがあれば追記します.