pari-gpで代数体の計算（2次体の続き）
はじめに
前回の記事pari-gpで代数体関連の計算をする の続きです．
代数体の計算と言いつつ，2次体，さらには虚2次体しか扱っていませんでした．
今回も，2次体の計算にまつわるいくつかの関数（qfb一族）を紹介し，同時にpari-gpでの簡単なプログラミングの例を挙げたいと思います．
アルゴリズムについて
pari-gpに実装されている，代数体関連のアルゴリズムは，H. Cohenさんの2冊の本 A Course in Computational Number Theory, Advanced Topics in Computational Number Theory, Springer-Verlag, に解説があります．
2次体のにまつわるいくつかの関数
前回の記事で紹介したbnfinit()をはじめとするbnf一族は，big number fieldつまり，代数体に関する多めのデータを扱う関数の一族です．（bnfは，Buchman's number fieldでもあります．代数体のイデアル類群，単数群などを計算するアルゴリズムを提案したBuchmanにちなむものです．Buchmanのsubexponential algorithmについては，Cohen上掲書をご覧ください）．
これに対して，計算が短時間で終わる整数底や判別式などのデータのみを計算する，nfinit()をはじめとするnf一族もあります．
nf一族，bnf一族とも，一般の代数体（有理数体$\mathbb{Q}$上の一変数既約多項式の根を$\mathbb{Q}$に添加した体）一般に適用されるアルゴリズムを実装したものです．一方で，2次体には，固有の事情を用いたより高効率のアルゴリズムがあります．
2次体に固有の一連の関数として，quad一族があります．また，2元2次形式に固有の一連の関数として，qfb一族（qfbはQuadratic Forms, Binaryの略だと思われますが，これは著者の推測です）も2次体と関連します．
計算結果はどのくらい正しいのか？
pari-gpでのイデアル類群，単数群などの計算は，どれもGRH (Generalized Riemann Hypothesis)を仮定してのものです．例えば，bnfinit()などの出力もそうです．この計算結果からGRHを取り除くための関数bnfcertify()が用意されています．
gp> nf=bnfinit(x^2+3299);
time = 3 ms.
gp> nf.no
%12 = 27
gp> nf.clgp
%13 = [27, [9, 3], [[3, 0; 0, 1], [11, 0; 0, 1]]]
gp> bnfcertify(nf)
time = 2 ms.
%14 = 1
bnfinit()の結果をbnfcertify()に与えて，1が返ってくればGRHなしで結果が正しいということです．そうでない場合，計算が終了しなかったり，メモリ不足などで異常終了した場合は，その検証ができなかったと言うことです．
このあたりの詳細については，users.pdfの3.13.7 Class group, units, and the GRHの節（pari-gp 2.11.1同梱，バージョンが変わるとセクション番号も変わる可能性があります）を参照して下さい．
類数とイデアル類群
前回の記事で，類群とイデアル類群を計算するのにbnfinit()を使いました．2次体ならquadclassunit()を使う方が効率がよいとされています．
gp> v=quadclassunit( -3299)
%42 = [27, [9, 3], [Qfb(3, 1, 275), Qfb(23, -17, 39)], 1]
(13:56) gp> v.no
%43 = 27
(13:56) gp> v.clgp
%44 = [27, [9, 3], [Qfb(3, 1, 275), Qfb(23, -17, 39)]]
quadclassunit()のパラメタは2次体の判別式です．この関数の返値は，やはりベクトル（数学的な意味ではなく，多くのプログラミング言語で固定長リストと呼ばれるデータ構造だと思って下さい）で，第1番目が類数，次がイデアル類群，次がイデアル類群の生成元（を2元2次形式）です．
.noで類数，.clgpでイデアル類群といった情報にアクセスできます．
上の例は，虚2次体のイデアル類群の3-rankが2になる最初の（判別式の絶対値が最小の）例です．（蛇足ながら，素数$p$と有限アーベル群$A$について，$A$の$p$-rankとは，$A/A^p$の$p$元体上の線型空間としての次元のことです）．
コンピュータを用いた探索・gpのif/for
数学の研究でコンピュータを使うときの典型的な使い方の一つとして，このように「ある条件を満たす対象を探す，数える」，といったことが挙げられます．上の例を探すには，次のようにしてforループで判別式をある範囲で動かしながら，ifや論理演算子を用いて条件を満たすものを使います．
gp> for(d=2,4000,if(isfundamental(-d), v=quadclassunit(-d); if(#v.cyc>1 && v.cyc[1]%3==0 && v.cyc[2]%3==0, print(-d, ", ", v))))
-3299, [27, [9, 3], [Qfb(3, 1, 275), Qfb(23, -17, 39)], 1]
-3896, [36, [12, 3], [Qfb(3, 2, 325), Qfb(31, 14, 33)], 1]
time = 177 ms.
for(X=a, b, seq)は，変数Xをaからbまで動かしながら，GPの文seqを評価し，最後の評価値を返します．if(a, {seq1}, {seq2})は文字通りifで，aが非零ならseq1が，aが零ならseq2が評価されます．
上の実行例だと，$d$を2から4000まで動かし，$-d$が基本判別式か否か（isfundamental()は引数が基本判別式なら1, そうでないなら0を返す関数）をみます．$-d$が基本判別式なら，quadclassunit()で虚2次体$\mathbb{Q}(\sqrt{-d})$のイデアル類群を計算しvに格納します．さらに，v.cycの長さが1より大か否かをみます．ベクトルvに対して#vでその長さが得られますので，#v.cycはイデアル類群のcyclic componentの長さが得られます．#v.cyc > 1ならイデアル類群は巡回群ではないということです．このとき，1番目と2番目の巡回成分の位数 v.cyc[1], v.cyc[2] が3で割り切れれば，イデアル類群の3-rankが2以上と言うことです．
gpの関数
上で書いた「イデアル類群の3-rankが2以上の虚2次体を列挙するfor文」を関数に仕立ててみましょう．
gp> enum_prank_two(a, b, p) = {my(rv); rv=List([]); for(d=a,b,if(isfundamental(-d), v=quadclassunit(-d); if(#v.cyc>1 && v.cyc[1]%p==0 && v.cyc[2]%p==0, listput(rv, [d,v])))); return(rv);}
gp> enum_prank_two(2,4000,3)
time = 156 ms.
%81 = List([[27, [9, 3], [Qfb(3, 1, 275), Qfb(23, -17, 39)], 1], [36, [12, 3], [Qfb(3, 2, 325), Qfb(31, 14, 33)], 1]])
GPの関数定義は，funcname(arg1, arg2, ...) = {body} の要領です．bodyの部分にはGPの文が来ます． enum_prank_twoが関数名，a, bは探索範囲で，$-b \le -d \le -a$の範囲の基本判別式$-d$で，$\mathbb{Q}(\sqrt{-d})$のイデアル類群の$p$-rankが2以上のものを探します．関数内で局所的な変数の定義は my() を使います（Perlを思い出す人もいるでしょう）．空ベクトル [] をリストにするのに List() を使っています．これで可変長のリストが得られ，条件にあう虚2次体のイデアル類群がみつかると，判別式とイデアル類群をペアにしたベクトルをrvに保存し（listput()），最後に返値として返します．
$p=5$の場合を試してみると，
gp> enum_prank_two(8000,16000,5)
time = 424 ms.
%89 = List([[11199, [100, [20, 5], [Qfb(2, 1, 1400), Qfb(40, 1, 70)], 1]], [12451, [25, [5, 5], [Qfb(5, 3, 623), Qfb(7, 3, 445)], 1]]])
となります（探索が8000からはじまっているのは，それより前に存在しないことを確認しているからです）．
類数が9で割り切れる虚2次体は沢山ありますが，ほとんどの虚2次体のイデアル類群は位数9の巡回群を部分群として含み，$3\times 3$型のアーベル群を部分群として含むものは稀です．実際に試してみると，
gp> c=0; for(d=2, 5000, if(isfundamental(-d) && qfbclassno(-d)%9==0,c=c+1)); c
time = 119 ms.
%1 = 114
gp> enum_prank_two(2, 5000, 3)
time = 275 ms.
%3 = List([[27, [9, 3], [Qfb(3, 1, 275), Qfb(23, -17, 39)], 1], [36, [12, 3], [Qfb(3, 2, 325), Qfb(31, 14, 33)], 1], [9, [3, 3], [Qfb(13, 9, 79), Qfb(17, 11, 61)], 1]])
gp> #%
%4 = 3
となります．最初の1行により，$-5000 \le -d \le 2$の範囲の基本判別式で，対応する虚2次体の類数が9で割り切れるものの個数が114個あることが分かります．一方，先ほど書いた enum_prank_two() を使うと，イデアル類群に$3\times 3$型のアーベル群を含む虚2次体は，同じ範囲に3個しかないことが分かります．
Cohen-Lenstra Heuristics/Arithmetic Statistics
位数9のアーベル群の同型類$A$は2つありますが，虚2次体のイデアル類群の部分群にあらわれる割合は $1/\lvert\text{Aut}(A)\rvert$ に比例するのではないか？という経験則から導かれる一連の予想があります．この経験則，またはそこから導かれる一連の予想を Cohen-Lenstra Heuristics と呼びますが，詳しいことはH. Cohen and H. Lenstra, Heuristics on Class Groups of Number Fields に譲ります．こういった，数論的に興味のある対象がどのように分布しているか，という問題は，近年Arithmetic Statistics の名前で大きく注目を集めています．
プログラムファイルの読み込み・データの読み書き
コマンドラインでpari-gpのプログラムを書くのは大変なので，ファイルとして書いておいて，それを読み込んで実行することを考えます．
例として，テキストエディタで（なんでもよいですがemacsがよいでしょう）ファイルに先ほどの enum_prank_two()を書いておきます．ファイル名は enum_prank_two.gp とでもしておきましょう．それをpari-gpに読み込むには
gp> \r enum_prank_two.gp
gp qfs = enum_prank_two(2, 4000, 3)
のようにできます．せっかく計算したものを，gpの終了と同時に失うのはもったいないので，次のようにファイルに書き出しておきます．
gp> write("qfs.txt", qfs)
gp>
これにより，計算結果がテキストファイルとして qfs.txt に書き出されます．いったんgpを終了しても，
gp> qfs = read("qfs.txt")
%10 = List([[27, [9, 3], [Qfb(3, 1, 275), Qfb(23, -17, 39)], 1], [36, [12, 3], [Qfb(3, 2, 325), Qfb(31, 14, 33)], 1], [9, [3, 3], [Qfb(13, 9, 79), Qfb(17, 11, 61)], 1]])
のように計算結果を復元することができます．
関数内で計算結果を書き出すことも，まったく同様です．
まとめ
今回も扱う対象は虚2次体に終始しましたが，pari-gpをスクリプト言語として扱う，はじめの部分に触れました．次回は，より代数的整数論っぽい話題を扱ってみたいと思います（が，どうなるか分かりません）．