pari-gpで代数体の計算2

# pari-gpで代数体の計算（2次体の続き）

## はじめに

 前回の記事[pari-gpで代数体関連の計算をする](https://qiita.com/iwaokimura/items/6877c8b01ace4e1d9dd4) の続きです．
 代数体の計算と言いつつ，2次体，さらには虚2次体しか扱っていませんでした．

今回も，2次体の計算にまつわるいくつかの関数（qfb一族）を紹介し，同時にpari-gpでの簡単なプログラミングの例を挙げたいと思います．

## アルゴリズムについて

pari-gpに実装されている，代数体関連のアルゴリズムは，H. Cohenさんの2冊の本 [A Course in Computational Number Theory](https://www.springer.com/gp/book/9783540556404), [Advanced Topics in Computational Number Theory](https://www.springer.com/gp/book/9780387987279), 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()``が用意されています．

```shell-session
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](http://pari.math.u-bordeaux.fr/pub/pari/manuals/2.11.1/users.pdf)の3.13.7 Class group, units, and the GRHの節（pari-gp 2.11.1同梱，バージョンが変わるとセクション番号も変わる可能性があります）を参照して下さい．

## 類数とイデアル類群

前回の記事で，類群とイデアル類群を計算するのに``bnfinit()``を使いました．2次体なら``quadclassunit()``を使う方が効率がよいとされています．

```shell-session
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``や論理演算子を用いて条件を満たすものを使います．
 
```shell-session
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文」を関数に仕立ててみましょう．

```shell-session
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$の場合を試してみると，

```shell-session
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$型のアーベル群を部分群として含むものは稀です．実際に試してみると，

```shell-session
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](http://hdl.handle.net/1887/2137) に譲ります．こういった，数論的に興味のある対象がどのように分布しているか，という問題は，近年[**Arithmetic Statistics**](http://swc.math.arizona.edu/aws/2014/index.html) の名前で大きく注目を集めています．

## プログラムファイルの読み込み・データの読み書き

コマンドラインでpari-gpのプログラムを書くのは大変なので，ファイルとして書いておいて，それを読み込んで実行することを考えます．

例として，テキストエディタで（なんでもよいですがemacsがよいでしょう）ファイルに先ほどの ``enum_prank_two()``を書いておきます．ファイル名は ``enum_prank_two.gp`` とでもしておきましょう．それをpari-gpに読み込むには

```shell-session
gp> \r enum_prank_two.gp
gp qfs = enum_prank_two(2, 4000, 3)
```

のようにできます．せっかく計算したものを，gpの終了と同時に失うのはもったいないので，次のようにファイルに書き出しておきます．

```shell-session
gp> write("qfs.txt", qfs)
gp>
```

これにより，計算結果がテキストファイルとして ``qfs.txt`` に書き出されます．いったんgpを終了しても，

```shell-session
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をスクリプト言語として扱う，はじめの部分に触れました．次回は，より代数的整数論っぽい話題を扱ってみたいと思います（が，どうなるか分かりません）．