2019-09-07

数理論理２：意味論

今回は
・言語，構造，割り当てに対して，その言語の項，論理式の値が定まること
そして
・言語に対して，その言語の valid な論理式の全体が定まること
を述べます．

以下，言語を固定します．

（定義２－１）構造は空でない集合Ｄ（domainと呼びます）と次のような写像Ｉ（解釈と呼びます）の対です．
　ｆがｎ項関数記号ならば
　　ｎ＞０ならば，Ｉ（ｆ）はＤ＾ｎからＤへの写像
　　ｎ＝０ならば，Ｉ（ｆ）∈Ｄ
　ｐがｎ項述語記号ならば
　　ｎ＞０ならば，Ｉ（ｐ）はＤ＾ｎから｛０，１｝への写像
　　ｎ＝０ならば，Ｉ（ｐ）∈｛０，１｝
また，Ｄの濃度＃Ｄを構造の濃度と呼びます．

なお，構造のことを解釈，或いは，モデルと呼ぶ流儀もあります．

次に，構造（Ｄ，Ｉ）を固定します．

（定義２－２）割り当てｓは変数記号全体からＤへの写像です．ｘ，ｙ，..．が変数記号，ｄ，ｅ，...∈Ｄならば，ｓのｘ，ｙ，..．の像をそれぞれｄ，ｅ，...に改めた割り当てをｓ［ｘ→ｄ，ｙ→ｅ，...］で表します．

（定義２－３）付値は（関数，述語，論理記号の出現の個数に関する帰納法により定義される）次のような写像｜｜です．
　ｓが割り当て，ｘが変数記号ならば，｜ｘ｜＝ｓ（ｘ）
　ｇがｎ項関数または述語記号，ｔ１，...，ｔｎが項ならば
　　ｎ＞０ならば，｜ｇ（ｔ１，...，ｔｎ）｜＝Ｉ（ｇ）（｜ｔ１｜，...，｜ｔｎ｜）
　　ｎ＝０ならば，｜ｇ｜＝Ｉ（ｇ）
　Ａ，Ｂが論理式ならば
　　｜￢Ａ｜＝１－｜Ａ｜
　　｜Ａ∧Ｂ｜＝｜Ａ｜×｜Ｂ｜
　　｜Ａ∨Ｂ｜＝１－（１－｜Ａ｜）×（１－｜Ｂ｜）
　　｜Ａ→Ｂ｜＝１－｜Ａ｜×（１－｜Ｂ｜）
　ｓが割り当て，ｘが変数記号，Ａが論理式ならば
　　｜∀ｘＡ｜＝｛ｓ［ｘ→ｄ］に対して定義された｜｜による｜Ａ｜：ｄ∈Ｄ｝の元の積
　　｜∃ｘＡ｜＝｛ｓ［ｘ→ｄ］に対して定義された｜｜による｜Ａ｜：ｄ∈Ｄ｝の元の和

この定義において，付値による各元の像を元の値と呼びます．各項の値はＤに属し，各論理式の値は｛０，１｝に属します．また，各割り当てに対して，付値は項と論理式全体の集合からＤ∪｛０，１｝への準同型写像になっています．

以下，言語のみを固定します．

（例２－１）Ａ，Ｂが論理式ならば，任意のＳ，ｓに対して
　　｜（（Ａ→Ｂ）→Ａ）→Ａ｜
　＝１－｜（Ａ→Ｂ）→Ａ｜×（１－｜Ａ｜）
　＝１－（１－｜Ａ→Ｂ｜×（１－｜Ａ｜））×（１－｜Ａ｜）
　＝１－（１－（１－｜Ａ｜×（１－｜Ｂ｜））×（１－｜Ａ｜））×（１－｜Ａ｜）
であり，｜Ａ｜∈｛０，１｝により，この値は１です．

（例２－２）ｘが変数記号，Ａが論理式ならば，任意のＳ，ｓに対して
　　｜￢（∀ｘＡ）｜
　＝１－｜∀ｘＡ｜
　＝１－（｛ｓ［ｘ→ｄ］に対して定義された｜｜による｜Ａ｜：ｄ∈Ｄ｝の元の積）
　＝｛１－（ｓ［ｘ→ｄ］に対して定義された｜｜による｜Ａ｜）：ｄ∈Ｄ｝の元の和
　＝｛ｓ［ｘ→ｄ］に対して定義された｜｜による｜￢Ａ｜：ｄ∈Ｄ｝の元の和
　＝｜∃ｘ（￢Ａ）｜
となります．

（定義２－４）Ｓが構造，ｓが割り当て，｜｜が付値，Ａ，Ｂが論理式，Σが論理式の（有限とは限らない）集合ならば
１．｜Ａ｜＝１を（Ｓ，ｓ）|=Ａと書き，「Ｓ，ｓはＡを充足」と読み，そうしたＳ，ｓが存在するようなＡを「充足可能」と呼びます．
２．任意のｓに対して（Ｓ，ｓ）|=ＡをＳ|=Ａと書き，「ＡはＳにおいてvalid」或いは「ＳはＡのモデル」と読みます．
３．任意のＢに対して，Ｂ∈ΣならばＳ|=ＢをＳ|=Σと書き，「ＳはΣのモデル」と読みます．
４．任意のＳに対して，Ｓ|=ΣならばＳ|=ＡをΣ|=Ａと書き，「ＡはΣの意味論的帰結」と読みます．
５．任意のＳ，ｓに対して，｜Ａ｜＝｜Ｂ｜をＡ≡Ｂと書き，「Ａ，Ｂは同値」と読みます．

（定義２－４）において

・　Ａが閉論理式ならば，｜Ａ｜は割り当てによらないので，｜Ａ｜＝１はＳ|=Ａであり（単に「ＳはＡを充足」と読みます），Ａが充足可能とＡがモデルを持つは等価です．
・　Σが有限集合ならば
　　　｛Ａ１，...，Ａｎ｝|=ＡはＡ１，...，Ａｎ|=Ａ
　　　｛｝|=Ａは|=Ａ
のように｛　｝を略して書くことがあり
　　　|=Ａを「Ａは valid」と読みます．
　　　Ａ|=Ｂは|=（Ａ→Ｂ）と等価であり，これをＡ⇒Ｂと書きます．
　　　Ａ|=ＢかつＢ|=Ａは|=（Ａ→Ｂ）∧（Ｂ→Ａ）と等価であり，これをＡ⇔Ｂと書きます．
　　　Ａ⇔ＢはＡ≡Ｂと等価です．
・　|=Ａは，任意のＳ，ｓに対して，（Ｓ，ｓ）|=Ａつまり｜Ａ｜＝１つまり｜￢Ａ｜＝０と書けるので，入力されたＡにおいて束縛されない全ての変数記号を定数記号とみなし，それがvalidであることを確かめよう，或いは，｜￢Ａ｜＝１となるＳ（これをＡのカウンターモデルと呼びます）を見付けてvalidでないことを確かめようとする https://www.umsu.de/trees/ のようなツールがあります．
・　論理式全体は，適当なＳ，ｓに充足される（充足可能といいます）論理式全体と，そうでない（充足不可能といいます）論理式全体に分割され，前者の部分集合として適当なＳにおいて valid な論理式全体があり，さらにその部分集合として valid な論理式全体があります．

（例２－１）により
　|=（（Ａ→Ｂ）→Ａ）→Ａ　つまり　（Ａ→Ｂ）→Ａ|=Ａ　つまり　（Ａ→Ｂ）→Ａ⇒Ａ
（例２－２）により
　￢（∀ｘＡ）≡∃ｘ（￢Ａ）　つまり　￢（∀ｘＡ）⇔∃ｘ（￢Ａ）
となります．

なお，（定義２－４）において，より強く

４．任意のＳ，ｓに対して，「任意のＢに対して，Ｂ∈Σならば（Ｓ，ｓ）|=Ｂ」ならば（Ｓ，ｓ）|=ＡをΣ|=Ａと書き，「ＡはΣの意味論的帰結」と読みます．

とする流儀もあります．この場合，Σ∪｛Ａ｝|=Ｂならば Σ|=Ａ→Ｂ（意味論的な演繹定理）がＡが閉論理式でなくても成り立つ（つまり，ｐが1項述語記号，ｘ，ｙが変数記号であっても｛ｐｘ｝|=ｐｙとならない）など，構文論との対応が失われますが，閉論理式に限定するなら差異はありません．

2019-09-07

nfsplitting(*,,2);

私が案出した Galois 群の計算方法を PARI の開発メンバーの方にお伝えしたところ，それに基づいた nfsplitting さらに galoisinit の改良版を作成してくださいました．試用方法は
https://pari.math.u-bordeaux.fr/cgi-bin/gitweb.cgi?p=pari.git;a=tree;h=refs/heads/bill-nfsplitting;hb=refs/heads/bill-nfsplitting
の snapshot のリンクをクリック，保存，展開して

./Configure;make all;./gp

で，PARI が起動します．以下，実行例です．

? #
   timer = 1 (on)
? galoisinit(nfsplitting(x^17-2)); \\ 従来の方法
  *** galoisinit: Warning: increasing stack size to 16000000.
  *** galoisinit: Warning: increasing stack size to 32000000.
  *** galoisinit: Warning: increasing stack size to 64000000.
time = 1h, 3min, 8,876 ms.
? galoissplittinginit(x^17-2); \\ 新作
  *** galoissplittinginit: Warning: increasing stack size to 16000000.
  *** galoissplittinginit: Warning: increasing stack size to 32000000.
  *** galoissplittinginit: Warning: increasing stack size to 64000000.
time = 20,935 ms.
? g=nfsplitting(f=x^17-2);nfisincl(f,g); \\ 従来の方法
  *** nfisincl: Warning: increasing stack size to 16000000.
  *** nfisincl: Warning: increasing stack size to 32000000.
time = 20,210 ms.
? nfsplitting(x^17-2,,1); \\ 処理の共通部分をまとめたもの．出力は分解体の定義多項式と入力の全ての根の分解体上での表示です．
  *** nfsplitting: Warning: increasing stack size to 16000000.
  *** nfsplitting: Warning: increasing stack size to 32000000.
time = 16,544 ms.
? nfsplitting(x^17-2,,2); \\ 出力は分解体の定義多項式とその Galois 群の写像表示です．
  *** nfsplitting: Warning: increasing stack size to 16000000.
  *** nfsplitting: Warning: increasing stack size to 32000000.
time = 16,846 ms.

2019-09-07

数理論理１：言語

ちょっとした経緯があり，数理論理（古典一階論理）の基本事項を少し書きます．

今回は
・言語に対して，その言語の項，論理式の全体が定まること
を述べます．

（定義１－１）言語は，次のような互いに区別できる記号の集合です．
　ｎ項関数記号　（各ｎは非負整数）
　ｎ項述語記号　（１個以上，各ｎは非負整数）
　変数記号　ｘ１　ｘ２　 …　（可算個）
　論理記号　￢　∨　∧　→　（結合子とも呼びます）　∀　∃　（量化子とも呼びます）
　補助記号　（　，　）

ここで，関数記号，述語記号は言語を特徴付けるもので
・それぞれの全体（と各記号に上記のｎを対応させる関数）の組をその言語のシグネチャーと呼びます．
なお
・関数記号，述語記号を非論理記号，その他（つまり，全ての言語に共通するもの）を論理記号と呼ぶ流儀，また，等号＝を論理記号に含める流儀もあります．

言語を固定すると，その言語の項，論理式が以下のように定まります．

（定義１－２）項は
　０項関数記号　（定数記号とも呼びます）
　変数記号
　適当なｎ（＞０）項関数記号ｆ，ｎ個の適当な項ｔ１，...，ｔｎ，補助記号によってｆ（ｔ１，...，ｔｎ）と表せるもの
のいずれかです．

（例１－１）ｘが変数記号，ｃが定数記号，ｆが２項関数記号ならば，ｆ（ｘ，ｃ）は項です．

なお，（定義１－２）では「逆にこれらはいずれも項です．」を略しており，正確には，現れる定数記号でない関数記号の個数に関する帰納法により

　Ｔ（０）＝変数記号の全体と定数記号の全体との和集合
　任意の非負の整数ｊに対して
　　Ｔ（ｊ＋１）＝∪＿｛ｎ≧１｝｛ｆ（ｔ１，...，ｔｎ）：ｆはｎ項関数記号，ｔ１，...，ｔｎはそれぞれＴ（ｍ＿１），...，Ｔ（ｍ＿ｎ）の元．ただし，ｍ＿１，...，ｍ＿ｎは和がｊである非負の整数．｝
と定め，項全体を
　∪＿｛ｊ≧０｝Ｔ（ｊ）
と定めます．

そして，ｆ（ｔ１，...，ｔｎ）は，項の個数に関する帰納法により

　ｆ（ｔ１，...，ｔ１）＝ｆ（ｔ１）
　任意の正の整数ｊに対して
　　ｆ（ｔ１，...，ｔ(ｊ＋１)）＝ｆ（ｔ１，...，ｔｊ）の最後の "）" を "，ｔ(ｊ＋１)）" に置き換えたもの
と定めます．

同様に（定義１－３），（定義２－３）は関数，述語，論理記号の出現の個数，（定義３－１）は推論規則を適用した回数に関する帰納法による定義です．

（定義１－３）論理式は
　０項述語記号　（命題記号とも呼びます）
　適当なｎ（＞０）項述語記号ｆ，ｎ個の適当な項ｔ１，...，ｔｎ，補助記号によってｐ（ｔ１，...，ｔｎ）と表せるもの
のどちらか（これらを原子論理式と呼びます），或いは，適当な変数記号ｘ，論理式Ａ，Ｂおよび論理記号，補助記号によって
　￢（Ａ）
　（Ａ）∧（Ｂ）
　（Ａ）∨（Ｂ）
　（Ａ）→（Ｂ）
　（∀ｘ）（Ａ）　（この∀ｘのｘ，および，Ａ内の各ｘをｘの（全称）束縛出現といいます）
　（∃ｘ）（Ａ）　（この∃ｘのｘ，および，Ａ内の各ｘをｘの（特称）束縛出現といいます）
のいずれかと表せるものです．

ここで
・論理式に現れる個々の変数記号をその変数記号の出現
・束縛出現でない変数記号の出現を自由出現
・各変数記号の出現がすべて束縛出現である論理式を閉（じた）論理式，そうでない論理式を開（いた）論理式
そして
・閉論理式以外の元をもたない集合を理論
と呼びます．

なお，Ａ，Ｂが論理式ならな
　（Ａ） $\leftrightarrow$ （Ｂ）は（（Ａ）→（Ｂ））∧（（Ｂ）→（Ａ））
の略記とします．

また，言語の定義において，論理記号は，例えば　￢　∧　∃　のみとして
　（Ａ）∨（Ｂ）は￢（（￢（Ａ））∧（￢（Ｂ）））
　（Ａ）→（Ｂ）は￢（Ａ∧（￢（Ｂ）））
　（∀ｘ）（Ａ）は￢（（∃ｘ）（￢（Ａ）））
それぞれの略記と定義する流儀もあります．

（例１－２）ｘ，ｙが変数記号，ｆが１項関数記号，ｐが１項述語記号，ｒが２項述語記号ならば
　ｒ（ｆ（ｘ），ｙ）
　（（∀ｘ）（ｐ（ｘ）））→（ｐ（ｘ））
はどちらも論理式，そして
　（∃ｘ）（（ｐ（ｘ））→（（∀ｘ）（ｐ（ｘ））））
は閉論理式です．

（例１－２）の論理式をそれぞれ
　ｒｆｘｙ
　（∀ｘｐｘ）→ｐｘ　或いは　（∀ｘ（ｐ（ｘ）））→ｐ（ｘ）
　∃ｘ（ｐｘ→∀ｘｐｘ）　或いは　∃ｘ（ｐ（ｘ）→∀ｘ（ｐ（ｘ）））
などのように記すこともあり，本稿もその作法に従います．

（定義１－４）Ａが論理式，ｘ，ｙ，..がＡに現れる束縛されていない全ての変数記号ならば，Ａの全称閉包は閉論理式∀ｘ∀ｙ...（Ａ）です．

論理式を"閉じる"方法としては，各変数の自由出現を各変数ごとに新たな定数記号に置き換えるというものもあり，本稿のメタ論理でも暗にそれを採用しています．

（参考）
https://www.encyclopediaofmath.org/index.php/Formalized_language

2019-07-18

代数体上の線型方程式系の求解

今回の冪根拡大の構成では，中間体上の線型方程式系（連立一次方程式）を解くことが一つの柱となっています（ http://ehito.hatenablog.com/entry/2019/04/11/193241 ）．それは例えば，定義多項式が

c5^4+c5^3+c5^2+c5+1,
e1^2+10,
e2^5+((-3750*c5^3)-3750*c5^2-1875)*e1-3125*c5^3+9375*c5^2+6250*c5+3125

である体 Q(c5,e1,e2) 上で

(((10*a4-2*a3+30*a2-70*a1+360)*e1-10*a4-20*a3-150*a2-100*a1-600)*c5^3
         +((10*a4-14*a3+30*a2-90*a1+360)*e1-30*a4-20*a3-50*a2-100*a1-1000)
          *c5^2+(((-16*a3)-160*a1)*e1-40*a4-200*a2-1600)*c5
         +(5*a4-8*a3+15*a2-80*a1+180)*e1-240*e2-20*a4-10*a3-100*a2+250*a1-800)
         /240,
        (((2*a4-20*a3-10*a2-140*a1+40)*e1-20*a4-60*a3-100*a2+100*a1-800)*c5^3
         +((4*a4-20*a3-20*a2-140*a1+80)*e1-20*a4-100*a2+200*a1-800)*c5^2
         +((6*a4-30*a2+120)*e1-60*a3+300*a1)*c5
         +(3*a4-40*a3-15*a2-120*a1+60)*e1-48*e2^2-60*a4-30*a3-100*a2+150*a1
         -2000)
         /48,
        -(((125*a4+10*a3+575*a2+750*a1+4900)*e1
         +150*a4-450*a3+250*a2-2750*a1+5000)
         *c5^3
         +((125*a4-30*a3+575*a2-250*a1+4900)*e1
          -250*a4-450*a3-1750*a2-2750*a1-11000)
          *c5^2+(((-20*a3)+500*a1)*e1-100*a4-1500*a2-6000)*c5
         +(100*a4-10*a3+100*a2+250*a1+3200)*e1+48*e2^3-50*a4-600*a3-750*a2
         -4000*a1-3000)
         /48,
        -(((450*a4+1000*a3+3750*a2+3000*a1+21000)*e1
         -1500*a3-5000*a2-7500*a1-10000)
         *c5^3
         +((400*a4+1000*a3+3000*a1+12000)*e1-5000*a3-5000*a2-10000*a1-10000)
          *c5^2+((850*a4+3750*a2+33000)*e1-6500*a3-17500*a1)*c5
         +(425*a4-1000*a3+1875*a2-1000*a1+16500)*e1+48*e2^4-2500*a4-3250*a3
         -12500*a2-8750*a1-100000)
         /48

の零点 [a4,a3,a2,a1] を（少なくとも一つ）求めるといった処理であり，ツールとしては，Nemo/Julia シリーズ AbstractAlgebra.jl の solve_left ( https://nemocas.github.io/AbstractAlgebra.jl/latest/matrix/#AbstractAlgebra.Generic.solve_left-Union{Tuple{T},%20Tuple{MatElem{T},MatElem{T}}}%20where%20T%3C:RingElem )

Welcome to Nemo version 0.14.1

Nemo comes with absolutely no warranty whatsoever


Welcome to 

    _    _           _        
   | |  | |         | |       
   | |__| | ___  ___| | _____ 
   |  __  |/ _ \/ __| |/ / _ \
   | |  | |  __/ (__|   <  __/
   |_|  |_|\___|\___|_|\_\___|
    
Version 0.6.2 ... 
 ... which comes with absolutely no warranty whatsoever
(c) 2015-2019 by Claus Fieker, Tommy Hofmann and Carlo Sircana


julia> MXl2=MXl[2];

julia> MX=matrix(MXl[4],MXl2,MXl2,MXl[1]);

julia> cx=matrix(MXl[4],1,MXl2,MXl[3]);

julia> @time xx=solve_left(MX,cx);
  1.407416 seconds (611.78 k allocations: 29.782 MiB, 1.75% gc time)

および，maxima の algebraic モードでの linsolve と fast_linsolve

(%i11) (map(lambda([s],tellrat(s[1])),DP),algebraic:on,linsolve(sAL1,AL0))$
Evaluation took 0.0040 seconds (0.0040 elapsed) using 1.843 MB.

(%i12) (map(lambda([s],tellrat(s[1])),DP),algebraic:on,fast_linsolve(sAL1,AL0))$
Assuming entries of type ANY-MACSYMA
Starting to solve.  There are 4 equations with 4 unknowns occurring.
The value of (SP-TYPE-OF-ENTRIES SP-MAT) is ANY-MACSYMA
The dimension of the solution space is 0
Evaluation took 0.0040 seconds (0.0040 elapsed) using 1.186 MB.

があります．

これらのうち solve_left は最近公開されたもので，高速 arithmetic で知られる Nemo シリーズゆえ相応のパフォーマンスが期待され，従来の hnf_with_transform と can_solve_left_reduced_triu の組み合わせよりは高速化していますが，サイズの大きな問題に対する処理速度は maxima の fast_linsolve が優位であり，例えば，47 分多項式を入力とする処理に現れる，文字列としてのサイズが 405801 バイト（上記の例は 1087 バイト）の方程式系の場合，

julia> @time xx=solve_left(MX,cx);
80.273864 seconds (65.62 M allocations: 31.598 GiB, 1.47% gc time)

に対し，無印の linsolve では

(%i14) (map(lambda([s],tellrat(s[1])),DP),algebraic:on,linsolve(sAL1,AL0))$
Evaluation took 458.1150 seconds (472.1600 elapsed) using 92401.623 MB.

そして，

(%i15) (map(lambda([s],tellrat(s[1])),DP),algebraic:on,fast_linsolve(sAL1,AL0))$
Assuming entries of type ANY-MACSYMA
Starting to solve.  There are 22 equations with 22 unknowns occurring.
The value of (SP-TYPE-OF-ENTRIES SP-MAT) is ANY-MACSYMA
The dimension of the solution space is 0
Evaluation took 54.8670 seconds (56.4100 elapsed) using 20716.948 MB.

となります．

2019-07-08

位数の取得

前回も触れたように，nfsplitting の処理時間は結果の多項式の次数を第二引数として与えることで一般に短縮されます．この性質を利用するため本プログラムでは，その次数，つまり，入力多項式 $f$ の $\mathbb{Q}$ 上の Galois 群 $G$ の位数を，予め用意したある範囲内の全ての cycle index polynomial (https://en.wikipedia.org/wiki/Cycle_index) の中から $G$ の cycle index polynomial $Z(G)$ を特定することで取得しています．

置換群の cycle index polynomial とは，その群に属する元の cycle type ( https://groupprops.subwiki.org/wiki/Cycle_type_of_a_permutation ) の頻度を変数名と次数，そして係数により，一つの多変数多項式で表したもので，群の位数の逆数は（1個しかない）単位元の cycle type の頻度（対応する単項式の係数）として現れます．参照する cycle index polynomial のデータは，次数が $37$ 以下の $f$ の Galois 群を全て含む，GAP の Transitive Groups Library transgrp ( https://www.gap-system.org/Packages/transgrp.html ) から抽出したもので，その制約から次数が $38$ 以上の $f$ についてはノーヒントで nfsplitting を実行しています．なお，同型でない群でも同じ cycle index polynomial を有するものが存在するため，群自体の特定には至らないことはよく目にする注意です．

実際の $f$ から $Z(G)$ を特定するには，前回も引用した「正の数 $M$ 以下で $f$ の判別式の因数ではない素数 $p$ を法とする $f$ の既約因子の次数型（どのような次数の因子がいくつずつあるか）の密度が， $M\to+\infty$ のとき， $G$ の同じ型の cycle type の頻度に収束する」という Chebotarev's density theorem を踏まえ， $p$ が値の小さいほうから 10^3 個目の素数となる，または，単位元の cycle type が 3 回現れるまで，既約分解を繰り返し，各次数型と密度から cycle index polynomial と同じ書式の多項式を作り，用意したデータのうち，この多項式との差の係数ベクトルの $l_2$ ノルムが最小となるものを $Z(G)$ としています．

以下，この処理を組み込んだ関数 nfsplitting2 と標準の nfsplitting の実行例です．

? for(i=1,30,nfsplitting2(x^i-2));
[0,1,1]
[0.03125000000,2,1]
[0.004801097394,6,1]
[0.001226218787,8,1]
[0.0008579881657,20,1]
[0.002494223903,12,1]
[0.0007927650227,42,1]
[0.002837500000,16,1]
[0.0004409526648,54,1]
[0.003434499314,40,1]
[7.133182736E-5,110,1]
[0.001070169231,48,1]
[0.0004472513280,156,1]
[0.002071109694,84,1]
[0.0006910009183,120,1]
[0.001790521626,64,1]
[0.001964397853,272,1]
[0.001106438745,108,1]
[0.0002551165675,342,1]
[0.001759224398,160,1]
[0.001657411051,252,1]
[0.0006244174148,220,1]
[0.001948888552,506,1]
[0.0008731821745,500,1]
[0.0008779652962,312,1]
[0.0003889834061,486,1]
[0.002439757082,336,1]
[0.0008391018743,812,1]
  *** nfsplitting: Warning: increasing stack size to 16000000.
time = 10,490 ms.

? for(i=1,30,nfsplitting(x^i-2));
  *** nfsplitting: Warning: increasing stack size to 16000000.
  *** nfsplitting: Warning: increasing stack size to 32000000.
  *** nfsplitting: Warning: increasing stack size to 64000000.
time = 3min, 24,135 ms.

GAP の組み込み関数 ProbabilityShapes は同じ原理で作動し，例えば

gap> ProbabilityShapes(x^16-x^15-x^14+4*x^13-3*x^12-x^11+5*x^10-7*x^9+2*x^8+4*x^7-5*x^6+5*x^5-x^4-2*x^3+2*x^2-2*x+1);time;
[ 1947 ]
111832
gap> TransitiveGroup(16,1947);
t16n1947
gap> Size(last);
7962624

のように位数を得ることも可能ですが，一般に

? galoisord(x^16-x^15-x^14+4*x^13-3*x^12-x^11+5*x^10-7*x^9+2*x^8+4*x^7-5*x^6+5*x^5-x^4-2*x^3+2*x^2-2*x+1)
time = 90 ms.
%2 = [0.006932446267,7962624,1]

のような高速処理は望めません．

2019-06-23

多項式の Galois 群計算

本プログラムの Galois 群計算の仕組みは http://ehito.hatenablog.com/entry/2019/02/24/200919 で述べましたが，今回は具体例を用いて解説します．この方法は，単に入力多項式の Galois 群 $G$ を特定するのではなく，本プログラムの目標である入力多項式の根（実質的には分解体の primitive element ）の冪根表示に必要な情報を取得する中で，根の置換と同型写像の像を構成するものであり，他の目的での利用価値は不明ながら，同じ結果を得るにあたり代数体上の因数分解を繰り返す方法よりも高速な処理が望めます（以下は基礎体が有理数体 $\mathbb{Q}$ の場合ですが，円分体でも方針は同じです）．

計算は次の５つのステップからなります．
（ステップ１）入力多項式 $f$ の $\mathbb{Q}$ 上の分解体の primitive element $\alpha$ の $\mathbb{Q}$ 上の最小多項式 $g$ を取得．
（ステップ２） $f$ の根を $\alpha$ に対応する変数 $A$ の $\mathbb{Q}$ 上の多項式で表したリスト $L$ を取得．
（ステップ３） $g$ の数値根（多倍長浮動小数点数）を $L$ の $A$ に代入し，根の置換としての $G$ の各元のリスト $G_1$ を取得．
（ステップ４）内積 $L\cdot K$ が $A$ となる整数のリスト $K$ を取得．
（ステップ５） $L\cdot K=A$ の $L$ に（ステップ３）で得た各置換を施し， $G$ の各元による $A$ の像のリスト $G_2$ を取得．

それでは， $f=x^{32}-2$ を例に各ステップを実行してみましょう．

（ステップ１）PARI/GP において

? d=poldegree(ga=subst(g=nfsplitting(f=x^32-2),x,A));g

と入力すると，約 10 秒（Ubuntu 14.04, Core i5-3210M 2.50GHz）で次のように出力されます（環境によっては default コマンドで parisizemax を増やす必要があるかも分かりません）．

%1 = x^256+9617286208*x^224+16939743610198607616*x^192+2420550419437751206410678272*x^160+5415888922248827561943116410019840*x^128+21997636389176028314448861495081828352*x^96+6678882252223940045087005013533786112*x^64+370736334465544641104295493632*x^32+16777216

なお，Chebotarev's density theorem ( https://en.wikipedia.org/wiki/Chebotarev%27s_density_theorem ) を用いて

? N=0;P=10^4;for(n=1,P,if(#factor(Mod(f,prime(n)))[,2]==poldegree(f),N=N+1));
time = 1,929 ms.
? nfsplitting(f,round(P/N));
time = 1,049 ms.

のように高速化できる場合もあります．
（ステップ２）次に

? L=nfisincl(f,ga);

と入力すると，やはり約 10 秒で結果が得られ，次数のみを拾うと

? apply(poldegree,L)
%3 = [241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241]

となっています．もし， $L$ を因数分解で得ようと

? nffactor(ga,f)

とすれば．．．time = 11min, 6,323 ms. を要します．
（ステップ３）続いて適当な精度を設定し， $g$ の数値根を求めます（精度不足ならば後の Map 生成時に一対一でない旨が指摘されます）．

? default(realprecision,max(200,floor(1.5*d)));ng=polroots(g);

それらを $L$ の $A$ に代入すると， $G$ の元の準同型性から，置換後の $f$ の数値根が得られるので，誤差を含む部分をカットします（精度不足の場合には，やはり Map 生成時に指摘されます）．

? nL=round(10^5*[subst(L,A,n)|n<-ng]);

この先頭のリストを基準に番号を割り振れば，根の置換群としての Galois 群 $G_1$ が得られます．

? M=Map(Mat([n1=nL[1]~,[1..#n1]~]));G1=apply(L->Vecsmall(apply(s->mapget(M,s),L)),nL);
time = 11 ms.

結果の中身を覗くと

? G1[256]
%10 = Vecsmall([32,3,2,28,29,7,6,24,25,11,10,20,21,15,14,16,17,19,18,12,13,23,22,8,9,27,26,4,5,31,30,1])
? G1[1]
%11 = Vecsmall([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32])

となっています．
（ステップ４） $L\cdot K=A$ を満たす $K$ の例を求めると

? K=matinverseimage(matconcat(vector(d,i,subst(L,A,i))~),[1..d]~)
%13 = [-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]~

となり，確かに

? L*K
%14 = A

となっています．
（ステップ５）（この例では不具合となりませんが， $G$ の位数の大きな問題では $L$ のサイズも大きくなるので）係数を逆に置換し， $A$ の像のリスト $G_2$ を計算します．

? G2=apply(p->L*p~,apply(p->apply(s->K[s],p),apply(p->Vec(Vecsmall(p)^(-1)),G1)));
time = 245 ms.

結果の中身を覗くと

? G2[1]
%16 = A
? G2[256]
%17 = -140056475292934525788582404664588050206670069515460087095141106743/35076811263842343718440191339547566716049104354229383482035976026041335790325202944*A^241-220739581799915126746789240596502331038977678547961/548075175997536620600627989680430729938267255534834116906812125406895871723831296*A^225-45346189340689199256448759270211222876248385931559781648808637129058905/1180878375432343917265021254361283554943748463312327749866549152506104759975936*A^209-2122915735604011586670819590599825988863612107748621607679979/548075175997536620600627989680430729938267255534834116906812125406895871723831296*A^193-37070637225162918656856440873437887413938164543428519572466260441003191806210881907/548075175997536620600627989680430729938267255534834116906812125406895871723831296*A^177-233704495019564084166805443788958628996577558892455075732292184313125/34254698499846038787539249355026920621141703470927132306675757837930991982739456*A^161-662136250029854892894980319171183536258815554834684143009196788040044982406610465898717709/68509396999692077575078498710053841242283406941854264613351515675861983965478912*A^145-8348613864254667748980481589124398754157263605209916044263166480616892223077/8563674624961509696884812338756730155285425867731783076668939459482747995684864*A^129-185188064704155447844565893173075231389171389924625227439295053172443812433789217646889315355429/8563674624961509696884812338756730155285425867731783076668939459482747995684864*A^113-583740749888698264024996895890648074047860624270764984232490179816787253679769389/267614832030047178027650385586147817352669558366618221145904358108835874865152*A^97-94021954878039241476007349158712536023216697732706026258402809031488357341671233165947721189077867/1070459328120188712110601542344591269410678233466472884583617432435343499460608*A^81-1185485609681620689614824041651980277946739037704685809931592407246459157786629070661/133807416015023589013825192793073908676334779183309110572952179054417937432576*A^65-3568346813747491847830869204262761414148872467371423771007232604647490123832556602374866005560817/133807416015023589013825192793073908676334779183309110572952179054417937432576*A^49-22495943806238745534925077213822835558283166332763585419821490732033000276871382989/8362963500938974313364074549567119292270923698956819410809511190901121089536*A^33-24759304369138635226898064742818791415078768965335903993868383116707204581031081330150511/16725927001877948626728149099134238584541847397913638821619022381802242179072*A^17-4141308750802347261186343968762150510786111593016503005587014426837918040267/2090740875234743578341018637391779823067730924739204852702377797725280272384*A

となっています．勿論， $G_2$ を因数分解で得ようと

? nffactor(ga,g)

とすれば．．．time = 1h, 38min, 46,267 ms. を要し，galoisinit や nfgaloisconj でも同様です．

なお，この $f$ の場合，冪根拡大を含む全体の処理時間は約 30 秒，そのうち約 25 秒が上記５つのステップに費やされます．一方，online の Magma V2.24-5 ( http://magma.maths.usyd.edu.au/calc/ ) の SolveByRadicals ( https://magma.maths.usyd.edu.au/magma/handbook/text/401 https://carma.newcastle.edu.au/tdlc/sins/190503_nicole_sutherland.pdf http://www.fachgruppe-computeralgebra.de/data/CA-Rundbrief/car62.pdf ) で

P<x> := PolynomialRing(IntegerRing());
SolveByRadicals(x^32 - 2);

のように処理すると

Runtime error: too many transformation, cannot find a primitive element

となります．

2019-05-22

円分体上の Galois 群

前回までの冪根拡大では，円分体を基礎体とする一方，中間体の生成は $\mathbb{Q}$ 上の Galois 群の組成列に従っており， $M_i(X,0)\in F_{i-1}[X]$ となる $i$ の発生原因となっていました．今回はこれを「まず円分体とその上の分解体の Galois 群を構成し，その組成列を用いて中間体を生成する手順」に，また，「最小多項式の文字列としてのサイズに応じて， $M_i(X,0)$ の展開方法を切り替える方式」（下の表示に出力があるものは各外部ツール，その他は組み込み関数を使用）に改めましたので，サンプルの実行結果を公開します．

表示は
・サンプル番号
・入力多項式
・ $\mathbb{Q}$ 上の分解体の primitive element の $\mathbb{Q}$ 上の最小多項式の次数
・同じく円分体上の最小多項式の（主係数についての）次数
・各中間体の拡大次数のリスト
・組成列までの処理時間（秒）
・冪根拡大の処理時間（秒）
のリストです．

展開の外部ツールとして Nemo を用いた場合．

[1,x^2-2,2,2,[2],0.152,0.002] 
[2,x^3-3*x-1,3,3,[3],0.137,0.003] 
[3,x^4-2,8,8,[2,2,2],0.184,0.006] 
[4,x^4+x^2-1,8,8,[2,2,2],0.163,0.007] 
[5,x^4-2*x^3+2*x^2+2,12,12,[3,2,2],0.181,0.03] 
[6,x^4+2*x^3+3*x^2+4*x+5,24,24,[2,3,2,2],0.275,0.125] 
[7,x^4+x+1,24,24,[2,3,2,2],0.234,0.089] 
[8,x^5-2,20,5,[5],0.209,0.002] 
[9,x^5-5*x+12,10,10,[2,5],0.205,0.022] 
[10,x^5+20*x+32,10,10,[2,5],0.192,0.027] 
[11,x^5+11*x+44,10,10,[2,5],0.166,0.026] 
[12,x^5+x^4-4*x^3-3*x^2+3*x+1,5,5,[5],0.172,0.009] 
[13,x^5+100*x^2+1000,20,5,[5],0.21,0.012] 
[14,x^6+x^3+1,6,3,[3],0.144,0.003] 
[15,x^6-2,12,6,[2,3],0.184,0.004] 
[16,x^5-5*x^3+5*x-5,20,10,[2,5],0.19,0.022] 
[17,x^6+x^3+7,18,9,[3,3],0.166,0.011] 
[18,x^6-3*x^4+1,24,24,[2,3,2,2],0.266,0.061] 
[19,x^6+x^4-9,24,24,[2,3,2,2],0.219,0.083] 
[20,x^6+x^4-8,48,48,[2,2,3,2,2],0.611,0.703] 
[21,x^7+7*x^3+7*x^2+7*x-1,14,7,[7],0.184,0.039] 
[22,x^7-14*x^5+56*x^3-56*x+22,21,7,[7],0.19,0.033] 
[23,x^7-2,42,7,[7],0.21,0.004] 
[24,x^8-2*x^6-x^4+7*x^2-5*x+1,16,16,[2,2,2,2],0.28,0.034] 
[25,x^9+x^8+3*x^6+3*x^3-3*x^2+5*x-1,18,18,[2,3,3],0.253,0.049] 
[26,x^12-3,24,12,[2,2,3],0.516,0.009] 
[27,x^7-14*x^5+56*x^3-56*x+22,21,7,[7],0.192,0.029] 
[28,
 x^15-x^14-14*x^13+13*x^12+78*x^11-66*x^10-220*x^9+165*x^8+330*x^7-210*x^6
     -252*x^5+126*x^4+84*x^3-28*x^2-8*x+1,15,15,[3,5],0.606,0.182]
  
<permutation group with 72 generators>
Welcome to Nemo version 0.13.5
  0.634780 seconds (798.89 k allocations: 42.184 MiB, 1.88% gc time)
  3.561398 seconds (3.36 M allocations: 1014.657 MiB, 2.55% gc time)
[29,x^6+x^4-x^2+5*x-5,72,72,[2,2,2,3,3],4.496,15.317] 
  *** nfisincl: Warning: increasing stack size to 16000000.
  *** write: Warning: increasing stack size to 16000000.
<permutation group with 13 generators>
Welcome to Nemo version 0.13.5
  0.729594 seconds (834.13 k allocations: 43.711 MiB, 8.42% gc time)
  0.795107 seconds (1.01 M allocations: 50.953 MiB, 1.18% gc time)
[30,2*(x+1)^13-x^13,156,13,[13],12.036,12.616] 
  *** nfisincl: Warning: increasing stack size to 16000000.
<permutation group with 49 generators>
Welcome to Nemo version 0.13.5
  0.632685 seconds (811.09 k allocations: 42.324 MiB, 9.44% gc time)
  0.994494 seconds (1.16 M allocations: 75.809 MiB, 1.63% gc time)
Welcome to Nemo version 0.13.5
  0.639144 seconds (811.09 k allocations: 42.326 MiB, 9.76% gc time)
  0.808816 seconds (1.01 M allocations: 51.011 MiB, 1.18% gc time)
[31,
 x^14+28*x^11+28*x^10-28*x^9+140*x^8+360*x^7+147*x^6+196*x^5+336*x^4-546*x^3
     -532*x^2+896*x+823,98,49,[7,7],5.772,33.549]
  
  *** nfisincl: Warning: increasing stack size to 16000000.
  *** write: Warning: increasing stack size to 16000000.
<permutation group with 75 generators>
Welcome to Nemo version 0.13.5
  0.635679 seconds (802.47 k allocations: 41.997 MiB, 1.84% gc time)
  3.301811 seconds (2.28 M allocations: 812.627 MiB, 2.59% gc time)
Welcome to Nemo version 0.13.5
  0.688764 seconds (802.47 k allocations: 41.997 MiB, 8.95% gc time)
  1.009562 seconds (1.17 M allocations: 82.609 MiB, 1.53% gc time)
Welcome to Nemo version 0.13.5
  0.690487 seconds (802.47 k allocations: 41.997 MiB, 8.90% gc time)
  0.809244 seconds (1.01 M allocations: 50.949 MiB, 1.17% gc time)
[32,
 x^15-470*x^13-305*x^12+71840*x^11+85357*x^10-4292700*x^9-3714805*x^8
     +119761820*x^7+25284495*x^6-1542190154*x^5+717324725*x^4+7178878600*x^3
     -5452953875*x^2-7998223215*x+4461221029,75,75,[3,5,5],15.793,132.7]
  
Evaluation took 180.3810 seconds (240.6020 elapsed) using 74460.010 MB.

展開の外部ツールとして Giac を用いた場合．

[1,x^2-2,2,2,[2],0.133,0.001] 
[2,x^3-3*x-1,3,3,[3],0.169,0.003] 
[3,x^4-2,8,8,[2,2,2],0.168,0.005] 
[4,x^4+x^2-1,8,8,[2,2,2],0.16,0.007] 
[5,x^4-2*x^3+2*x^2+2,12,12,[3,2,2],0.174,0.029] 
[6,x^4+2*x^3+3*x^2+4*x+5,24,24,[2,3,2,2],0.273,0.126] 
[7,x^4+x+1,24,24,[2,3,2,2],0.224,0.087] 
[8,x^5-2,20,5,[5],0.182,0.002] 
[9,x^5-5*x+12,10,10,[2,5],0.18,0.023] 
[10,x^5+20*x+32,10,10,[2,5],0.163,0.032] 
[11,x^5+11*x+44,10,10,[2,5],0.146,0.025] 
[12,x^5+x^4-4*x^3-3*x^2+3*x+1,5,5,[5],0.141,0.008] 
[13,x^5+100*x^2+1000,20,5,[5],0.161,0.012] 
[14,x^6+x^3+1,6,3,[3],0.139,0.002] 
[15,x^6-2,12,6,[2,3],0.167,0.004] 
[16,x^5-5*x^3+5*x-5,20,10,[2,5],0.171,0.023] 
[17,x^6+x^3+7,18,9,[3,3],0.162,0.009] 
[18,x^6-3*x^4+1,24,24,[2,3,2,2],0.243,0.061] 
[19,x^6+x^4-9,24,24,[2,3,2,2],0.218,0.083] 
[20,x^6+x^4-8,48,48,[2,2,3,2,2],0.612,0.702] 
[21,x^7+7*x^3+7*x^2+7*x-1,14,7,[7],0.17,0.04] 
[22,x^7-14*x^5+56*x^3-56*x+22,21,7,[7],0.175,0.028] 
[23,x^7-2,42,7,[7],0.185,0.004] 
[24,x^8-2*x^6-x^4+7*x^2-5*x+1,16,16,[2,2,2,2],0.274,0.033] 
[25,x^9+x^8+3*x^6+3*x^3-3*x^2+5*x-1,18,18,[2,3,3],0.241,0.048] 
[26,x^12-3,24,12,[2,2,3],0.52,0.009] 
[27,x^7-14*x^5+56*x^3-56*x+22,21,7,[7],0.177,0.029] 
[28,
 x^15-x^14-14*x^13+13*x^12+78*x^11-66*x^10-220*x^9+165*x^8+330*x^7-210*x^6
     -252*x^5+126*x^4+84*x^3-28*x^2-8*x+1,15,15,[3,5],0.605,0.189]
  
<permutation group with 72 generators>
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
Added 0 synonyms

Evaluation time: 8.8
"Done","Done","Done","Done"
// Time 8.8
// Total time 8.8
[29,x^6+x^4-x^2+5*x-5,72,72,[2,2,2,3,3],4.448,15.771] 
  *** nfisincl: Warning: increasing stack size to 16000000.
  *** write: Warning: increasing stack size to 16000000.
<permutation group with 13 generators>
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
Added 0 synonyms
"Done","Done","Done","Done"
// Time 0
// Total time 0
[30,2*(x+1)^13-x^13,156,13,[13],12.067,8.419] 
  *** nfisincl: Warning: increasing stack size to 16000000.
<permutation group with 49 generators>
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
Added 0 synonyms

Evaluation time: 0.43
"Done","Done","Done","Done"
// Time 0.43
// Total time 0.43
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
Added 0 synonyms
"Done","Done","Done","Done"
// Time 0
// Total time 0
[31,
 x^14+28*x^11+28*x^10-28*x^9+140*x^8+360*x^7+147*x^6+196*x^5+336*x^4-546*x^3
     -532*x^2+896*x+823,98,49,[7,7],5.794,25.028]
  
  *** nfisincl: Warning: increasing stack size to 16000000.
  *** write: Warning: increasing stack size to 16000000.
<permutation group with 75 generators>
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
Added 0 synonyms

Evaluation time: 8.13
"Done","Done","Done","Done"
// Time 8.13
// Total time 8.13
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
Added 0 synonyms
"Done","Done","Done","Done"
// Time 0.35
// Total time 0.35
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
Added 0 synonyms
"Done","Done","Done","Done"
// Time 0.01
// Total time 0.01
[32,
 x^15-470*x^13-305*x^12+71840*x^11+85357*x^10-4292700*x^9-3714805*x^8
     +119761820*x^7+25284495*x^6-1542190154*x^5+717324725*x^4+7178878600*x^3
     -5452953875*x^2-7998223215*x+4461221029,75,75,[3,5,5],15.858,123.12]
  
Evaluation took 179.0620 seconds (218.4700 elapsed) using 74669.995 MB.