Big

明けましておめでとうございます。今年も宜しくお願い致します。

前回、素数を求めるC語をGo語に変換した事から、素数の特集をやってたニュートンを 借りてきて読んでる。

素数の定義は2以上の整数のうち1と自分自身でしか割る事が出来ない数って定義されてる。 ふーん、1はどうよ。上記を満たしていそうだぞ。なんで2からと断ってるの?

素数以外の数は、素数をかけて合成出来る。この合成数を作る素数の組み合わせは一通り しかない。しかし、1を素数と認定するとこれが破綻する。だから、1を除外するとな。

4以上の全ての偶数を2つの素数の和で表す事が出来るというのが、ゴールドバッハ予想と いう難問らしい。まだ、証明されていない。コンピュータの馬鹿力を借りて、400兆以下の 偶数について正しい事が確認されてるそうだ。

  16 =   5 +  11
 100 =  53 +  47
1000 = 821 + 179

素数を生成する関数は有るか? オイラーさんは、

f(x) = x^2 - x + 41

ってのを考えたけど、x=41 の所で破綻した。それじゃ、ある数pが素数かどうか判定する 方法は有るか? ウィルソンの定理というのが有るそうな。

! (p - 1) mod p == p - 1

なら、pは、素数。階乗の計算でっせ!!!

ガウスさんは、ある数Xまでの素数の数を予言した、素数定理って名前が付いているぐらい だから、証明されてるのかな?

π(x) = x / ln x

lnって、自然対数の事ね。

あと、世紀の大問題としてリーマン予想ってのが有るそうな。全ての整数で出来ている ゼータ関数と全ての素数で出来ているオイラー積が等しいってのの証明。素数に人生を 捧げて一生を終わる人は居るのだろうか?

math/big

去年はGo言語で年を締めたけど、今年もGoで行きます。新年なんで、Bigに行きます。 math/bigです。パッケージのカタログを眺めていて、たまたま見つけたやつ。 数字の大きいのを扱えるやつらしい。godoc math/bigすると、冒頭に

    Package big implements multi-precision arithmetic (big numbers). The
    following numeric types are supported:

        - Int   signed integers
        - Rat   rational numbers

ふむふむ、符号付きの整数と、有理数が扱えるとな。しかも精度は無限大。これはもう、 goでlispを組み立ててちょってためのパッケージか。それともLisp業界に喧嘩売るための 兵器かな?ああ、ruby業界も忘れてもらっちゃ困ると文句が出そうですけど。Java業界も、 俺んとこも有るぞと名乗りを挙げてきたりして。結局、近代的な言語は、巨大数とgcは必須 と言う事だな。

所で、このパッケージはどうやって使うの? そんな疑問は、上に続いて書いてありました。

    Methods are typically of the form:

        func (z *Int) Op(x, y *Int) *Int        (similar for *Rat)

    and implement operations z = x Op y with the result as receiver; if it
    is one of the operands it may be overwritten (and its memory reused). To
    enable chaining of operations, the result is also returned. Methods
    returning a result other than *Int or *Rat take one of the operands as
    the receiver.

math/bigのサンプル

ネットを探してみたけど、例が出てこなかったので、自分で例してみる。巨大数のサンプルと 言うと、factって決まってますな。再帰で書くのが普通だけど、ちょいと癖のあるインターフェースに なってるので、ループで回します。

package main

import (
        "fmt"
        "math/big"
        "os"
)

func main() {
        n := new(big.Int)
        _, err := fmt.Sscan(os.Args[1], n)
        if err != nil {
                panic(err)
        }
        one := big.NewInt(1) // one = 1
        r := new(big.Int)    // var r big.Int
        r.Set(n)             // r = n
        for {
                n.Sub(n, one)        // n = n - 1
                r.Mul(r, n)          // r *= n
                if n.Cmp(one) == 0 { // if n == 1
                        break
                }
        }
        fmt.Println(r)
}

実行例

sakae@uB:~/godev/src/t$ go run big.go 100
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

検算はいつものSchemeで行いました。合っててよかったね。ここまでは良いんですが、この パッケージはどう実装されてるの?

ソース読み

折角なんで、ソース読み(と言う観光、、、ですから、眺めるのが主目的)してみます。 もしemacsを使っているなら、上記のソースを開いておいて、例えばMulの所にカーソルを 置いて、M-. するだけ。(戻りはM-,) 使われているパッケージ(C語で言うとlib)の ソースが標準添付でさっと見れるのはありがたい。Cじゃこうはいかないものね。

// Mul sets z to the product x*y and returns z.
func (z *Int) Mul(x, y *Int) *Int {
    // x * y == x * y
    // x * (-y) == -(x * y)
    // (-x) * y == -(x * y)
    // (-x) * (-y) == x * y
    z.abs = z.abs.mul(x.abs, y.abs)
    z.neg = len(z.abs) > 0 && x.neg != y.neg // 0 has no sign
    return z
}

int.goに画面が移って、上記が表示されました。絶対値を先に計算して、後で符号を 調整してるみたい。本当の乗算は、パッケージ内だけでアクセス出来るmulに任せている んだな。(大文字で始まる関数は、公開。それ以外は非公開ってルールによる)

引き続いて、mulを見ると

func (z nat) mul(x, y nat) nat {
    m := len(x)
    n := len(y)

    switch {
    case m < n:
        return z.mul(y, x)
    case m == 0 || n == 0:
        return z.make(0)
    case n == 1:
        return z.mulAddWW(x, y[0], 0)
    }
    // m >= n > 1
     :

今度はnat.goなんて所に舞台が移り、数値の大きさで場合分けして再帰呼び出しとかで 単純にしてってる。mulAddWWを追って行くと、arith_decl.go で、

func mulAddVWW(z, x []Word, y, r Word) (c Word)

に到達したぞ。所で、上で出てきたnatって何よ?

type Int struct {
    neg bool // sign
    abs nat  // absolute value of the integer
}
type nat []Word
type Word uintptr

こんな風に、データ構造が定義されてたぞ。そんじゃ、最後に出てきたmulAddVWWはどうなってるの? さすがに、ここまで来ると、emacsにもリミッターがかかってるみたいで、見ちゃ駄目って 事みたい。見ると目が潰れるから。

math/bigのソースはどうなってるかツリーを辿ると

sakae@uB:/usr/lib/go/src/pkg/math/big$ ls
arith_386.s    arith.go           gcd_test.go      nat.go
arith_amd64.s  arith_test.go      hilbert_test.go  nat_test.go
arith_arm.s    calibrate_test.go  int.go           rat.go
arith_decl.go  example_test.go    int_test.go      rat_test.go

なんやらSが居ますなあ。おいらの石は糞386なので、arith_386.sを開くと

// func mulAddVWW(z, x []Word, y, r Word) (c Word)
TEXT ・mulAddVWW(SB),NOSPLIT,$0
        MOVL z+0(FP), DI
        MOVL x+12(FP), SI
        MOVL y+24(FP), BP
        MOVL r+28(FP), CX       // c = r
        MOVL z_len+4(FP), BX
        LEAL (DI)(BX*4), DI
        LEAL (SI)(BX*4), SI
        NEGL BX                 // i = -n
        JMP E5

L5:     MOVL (SI)(BX*4), AX
        MULL BP
          :

落ちるとこまで、落ちたな。結局、ユーザーに対応するのは、int.goでその下請けがnat.go 最下位は、arith*.sと言う構成なんだな。じゃ、これをどうやって組み立ててる?

sakae@uB:~/godev$ go install -x math/big
WORK=/tmp/go-build391788799
mkdir -p $WORK/math/big/_obj/
mkdir -p $WORK/math/
cd /usr/lib/go/src/pkg/math/big
/usr/lib/go/pkg/tool/linux_386/8g -o $WORK/math/big/_obj/_go_.8 -p math/big -D _/usr/lib/go/src/pkg/math/big -I $WORK ./arith.go ./arith_decl.go ./int.go ./nat.go ./rat.go
/usr/lib/go/pkg/tool/linux_386/8a -I $WORK/math/big/_obj/ -o $WORK/math/big/_obj/arith_386.8 -D GOOS_linux -D GOARCH_386 ./arith_386.s
/usr/lib/go/pkg/tool/linux_386/pack grcP $WORK $WORK/math/big.a $WORK/math/big/_obj/_go_.8 $WORK/math/big/_obj/arith_386.8
mkdir -p /usr/lib/go/pkg/linux_386/math/
cp $WORK/math/big.a /usr/lib/go/pkg/linux_386/math/big.a
go install math/big: open /usr/lib/go/pkg/linux_386/math/big.a: permission denied

なる程。goファイル群をコンパイル、sファイルをアセンブル。それらをパックしたのが、 パッケージの正体なのね。そして、それを所定の場所に置くとな。(権利が無いので最後の 所は、エラーだけど)

ウブに入ってるgoは、ちと古いので、言われた事を素直に実行するけど、1.3系だと利口に なってて、ソース群とパッケージのバイナリーファイルの作成日時を確認し、バイナリーが 新しければ何もしない、make風の挙動になってる。だから、コンパイルを発動させるには、 touch int.goとかして、goをだませば良い。

多倍長整数の演算

上記でざっとmath/bigを見たけれど、ちょっと消化不良。何か資料は無いものかとぐぐる先生に 聞いてみた所、アルゴリズムについて紹介多倍長整数の演算を教えてもらった。

加減乗除は、大きな桁数をまとめて、それを1桁として扱い、筆算の要領で演算してるのね。 納得しました。math/bigの表現だと、Vってのが1桁で、それが並んだのがWと表現してるんだな。 そして、符号の事は後回しにして、絶対値で計算してるとな。

多倍長整数を高速で行うためフーリエ変換法が有るって聞いていたけど、 高速フーリエ変換を利用した乗算処理の高速化 を見るまで、原理が分からなかった。理解に一歩近づいたな。それにしても、いろいろと 考える人がいるもんだ。

ああ、何年かぶりに取り出した、アルゴリズム辞典にも例が載ってたぞ。

gcd

再び math/bigに戻って上のソース群を見てたら、gcdなんて語句が目に留まった。gcdって 何だっけ? SICPにも出てきたような。。。あやつの1章は数学だからなあ。そうか、数学の 復習が必要か。以前見つけておいたサイト 高校数学+α :基礎と論理の物語を再訪。 第1章に出てた。

最大公約数(greatest common divisor)なんだな。これって、有理数(分数)を扱う時の約分に 使うのかな。まてまて、他にも使い道があるぞ。素数関係ね。高校数学でさりげなく素数を 最大限に利用したRSA暗号が説明されてた。これも時代の流れ?

素数の判定

math/bigをつらつら見てたら

func (x *Int) ProbablyPrime(n int) bool
    ProbablyPrime performs n Miller-Rabin tests to check whether x is prime.
    If it returns true, x is prime with probability 1 - 1/4^n. If it returns
    false, x is not prime.

多分素数だろうって判定してくれる関数も提供されてた。どんな事やってるか軽く見てみると、 nat.goに本体が鎮座してて、軽く素数か確認。後は、ごにょごにょ。

こんなのを見つけました。 Pythonを使って高速素数判定をしてみる これですっきりかな。

そして、こういう突っ込みも出てたぞ。 RSA鍵の生成時に確率的素数判定法を使って問題ないのか

math/bigの使い道

大体、どこで使われているか、想像が付くけど、地道に探ってみる。 cryptoって暗号パッケージで多用されてた。それから、encoding/asn1の中ね。大きな桁の 素数が必要な所。素数かどうか判定してるだろうから、

sakae@uB:/usr/lib/go/src/pkg/crypto$ grep ProbablyPrime */*.go
dsa/dsa.go:             if !q.ProbablyPrime(numMRTests) {
dsa/dsa.go:                     if !p.ProbablyPrime(numMRTests) {
rand/util.go:           if p.ProbablyPrime(20) && p.BitLen() == bits {
rsa/rsa.go:     // ProbablyPrime are deterministic, given the candidate number, it's
rsa/rsa.go:             if !prime.ProbablyPrime(20) {

この中で、素数を作っているんだな。で、現実的にはhttpsの通信に使うとな。sslがよく 使われているけど、Goが用意してるのは、tlsってやつ。ここでも、crypto/rsaをimport してるから、桁数の多い演算を使ってるって事になる。

以上から、math/bigは、近代のWebをサポートするのに必須のパッケージであり、決して Lispに対抗する為にサポートしてる訳では無い事が分かる。と、まあ、月並みな結論でした。

ssl

ここからはGoと離れ、巷でよく使われている(資料が多い)SSLをちと勉強してみる。

まずは、opensslの基本的な使い方って事で、 暗号化・復号化・公開鍵などを扱うツール

そして、素のツールを使って opensslでRSA暗号と遊ぶ 事が出来ますって。

最後は、Webでそれらがどう利用されてるか。 HTTPSコネクションの最初の数ミリ秒 の挙動の説明。

折角なんでsslとかのソースを観賞してみる。ソースの館、FreeBSDにスイッチします。

[sakae@fb10 /usr/src/crypto]$ ls
README  heimdal openssh openssl
[sakae@fb10 /usr/src/crypto/openssl/crypto/bn]$ ls
Makefile        bn_depr.c       bn_lib.c        bn_print.c      bntest.c
asm             bn_div.c        bn_mod.c        bn_rand.c       divtest.c
bn.h            bn_err.c        bn_mont.c       bn_recp.c       exp.c
bn.mul          bn_exp.c        bn_mpi.c        bn_shift.c      expspeed.c
bn_add.c        bn_exp2.c       bn_mul.c        bn_sqr.c        exptest.c
bn_asm.c        bn_gcd.c        bn_nist.c       bn_sqrt.c       todo
bn_blind.c      bn_gf2m.c       bn_prime.c      bn_word.c
bn_const.c      bn_kron.c       bn_prime.h      bn_x931p.c
bn_ctx.c        bn_lcl.h        bn_prime.pl     bnspeed.c

opensslの中の暗号の部のbnってのは、BigNumberの略なんでしょうかね。どんな構成に なってるか俯瞰してみると、

[sakae@fb10 /usr/src/crypto/openssl/crypto/bn]$ ls
Makefile        bn_depr.c       bn_lib.c        bn_print.c      bntest.c
asm             bn_div.c        bn_mod.c        bn_rand.c       divtest.c
bn.h            bn_err.c        bn_mont.c       bn_recp.c       exp.c
bn.mul          bn_exp.c        bn_mpi.c        bn_shift.c      expspeed.c
bn_add.c        bn_exp2.c       bn_mul.c        bn_sqr.c        exptest.c
bn_asm.c        bn_gcd.c        bn_nist.c       bn_sqrt.c       todo
bn_blind.c      bn_gf2m.c       bn_prime.c      bn_word.c
bn_const.c      bn_kron.c       bn_prime.h      bn_x931p.c
bn_ctx.c        bn_lcl.h        bn_prime.pl     bnspeed.c

大きな桁数の素数を求めるやつが有るなあ。bn_sqrt.cは素数の検索範囲の決定に使われる のかな。ぱある版の素数も求めるやつも混じってる。中を覗いてみると、2048個の素数を 発生してるぞ。ああ、これは、bn_prime.h内の素数テーブルを作成するのに使ってるのか。

詳しくは、 素数生成アルゴリズムの調査・開発 調査報告書 を見ろとな。ああ、思えば、遠くに来たもんだ。

おまけ

[sakae@fb10 ~]$ openssl prime 31
1F is prime
[sakae@fb10 ~]$ openssl prime 331
14B is prime
[sakae@fb10 ~]$ openssl prime 3331
D03 is prime
[sakae@fb10 ~]$ openssl prime 33331
8233 is prime
[sakae@fb10 ~]$ openssl prime 333331
51613 is prime
[sakae@fb10 ~]$ openssl prime 3333331
32DCD3 is prime
[sakae@fb10 ~]$ openssl prime 33333331
1FCA053 is prime
[sakae@fb10 ~]$ openssl prime 333333331
13DE4353 is not prime