Pony 試乗記(９) 多倍長演算ライブラリ作成

はじめに

Pony では倍精度浮動小数点数の他に 128bit 整数まで使えるので一般的な用途では困ることはないと思いますが、もっと高精度で計算したい、という場合には多倍長精度演算ライブラリが使えると便利です。

多倍長演算ライブラリの GMP (Gnu Multiple Precision Arithmetic Library) を Pony から扱えるようにしてみました。

動作確認のために、円周率を計算してみました。

方針

チュートリアルにもあるように、Pony では +, -, *, /,... のオペレータをメソッドとして定義できます。

GMP の mpf_t をラップしたクラスを作って、そのクラスに add, sub, mul, div, ... を実装すれば、通常のプリミティブな整数や浮動小数点数と同様に数式の形で記述できて便利です。

というわけで、概略次のようなクラスを作ることにしました:

class Mpf
  """
  Mpf クラスの概略
  """
  let _f: Pointer[None]  // mpf_t のアドレス

  new iso create() =>
    _f = @pony_mpf_init()  // mpf_t 領域を確保すると同時に、mpf_init を呼び出す

  fun _final() =>
    @pony_mpf_clear(_f)    // mpf_clear を呼んだうえで、mpf_t 領域を解放する

  (略)

  fun box add(other: Mpf val): Mpf val =>
    var r: Mpf iso = recover iso Mpf end  // 戻り値の mpf_t を作る
    @__gmpf_add(r._f, _f, other._f)       // mpf_add を呼び出す
    consume r                             // mpf_add の結果を返す

  fun box sub(other: Mpf val): Mpf val =>
    var r: Mpf iso = recover iso Mpf end  // 戻り値の mpf_t を作る
    @__gmpf_add(r._f, _f, other._f)       // mpf_add を呼び出す
    consume r                             // mpf_add の結果を返す
  (略)

これを使って、以下のように円周率を計算することができます。

  fun calc_pi(loop: I32 = 100): Mpf val =>
    """
    Gauss–Legendre algorithm
    """
    let f0: Mpf val = Mpf
    let f1: Mpf val = Mpf.from_i64(1)
    let f2: Mpf val = Mpf.from_i64(2)
    let f4: Mpf val = Mpf.from_i64(4)
    var a: Mpf val = f1
    var ba: Mpf val = f0
    var b: Mpf val = f1 / f2.sqrt()
    var t: Mpf val = f1 / f4
    var p: Mpf val = f1
    var k: I32 = 0
    while k < loop do
      ba = a
      a = (ba + b) / f2
      b = (ba * b).sqrt()
      t = t - ((p * (ba - a)) * (ba - a))
      p = f2 * p
      k = k + 1
    end
    ((a + b) * (a + b)) / (t * f4)

このコードは円周率を求めるプログラムを書くの記事を見て書きました。見比べて頂ければそのままなのがわかると思います。

[追記] 変数の Reference Capability はデフォルトで ref 扱いとなるので val より ref の方が書きやすいかと思い当初は以下のように定義していたのですが、ref だとアクターをまたがって値を渡すときなどで面倒なので、val に変更しました。

  new create() =>
    _f = @pony_mpf_init()

  fun box sub(other: Mpf val): Mpf val =>
    let r: Mpf = Mpf
    @pony_mpf_sub(r._f, _f, other._f)
    r

グルーコード

Pony からは引数にポインタを与えることはできますが実体を与えることはできません（そのうちできるようになるかも、とチュートリアルには書いてあります）。

一方で GMP の提供している API は、mpf_t 実体の操作として定義されています。

しかし、よくよく GMP のヘッダファイルを見ると共有ライブラリ上にある関数はポインタ渡しの形で記述されています。つまり、GMP のドキュメントにある

Function: void mpf_init (mpf_t x)

と呼び出したつもりでも、実際にはマクロで変換されて

typedef __mpf_struct *mpf_ptr;

void mpf_init (mpf_ptr);

が呼び出されることになっています。

なるべく GMP との手離れをよくしたかった（GMP の内部構造になるべく依存しない形にしたかった）ので、当初は C で公開 API を扱うグルーコードを書いて Pony からはそれを呼び出す形にしていたのですが、グルーコードを使うと環境依存でビルドできない可能性も出てきます。

mpf_t の構造の環境依存や、GMP のバージョンが上がって mpf_ptr を扱っている関数が変更される可能性と、グルーコードのビルドの手間とを考えた結果、グルーコードをなくして Pony から直接 mpf_ptr を扱う共有ライブラリ内の関数を呼び出す形にしました。

グルーコードをなくすにあたって、mpf_t の領域確保を Pony で行う必要が出てきました。このため、 (LGPL に汚染されないのか若干不安を感じつつも）mpf_t を struct として Pony 上で定義しました。

Corral 用パッケージ作成

折角なので、ある程度は人様に使っていただけるかたちにしたいと思いました。

以下の手順でパッケージにしてみました:

Library Scaffolding Generator を使ってパッケージの雛形を作成
作った雛形に、作成した GMP ラッパーソースを配置
corral.json にパッケージ情報を記述すると共に、post_fetch_or_update で make all を実行するように記述（共有ライブラリを自動的にビルドするため）

こうして作ったものは https://github.com/tadashi9e/gmp4pony.git にあります。

GMP ライブラリを利用するプロジェクトのディレクトリを作る。
プロジェクトディレクトリ上で corral init を実行する。
corral add https://github.com/tadashi9e/gmp4pony.git --version 0.0.1 を実行して取得先を設定する。
corral update を行う。このときに gmp4pony が取得されて共有ライブラリのビルドやテストケースの実行が行われるはず。
プロジェクト内の Pony ソースコードに use "gmp" を記述すれば Mpf クラスや Mpz クラスが利用可能になるはず。
corral run -- ponyc でプロジェクトをコンパイルする

という手順で利用可能になるはずです。

手元にある Ubuntu マシンでしか動作確認していないので Windows とかではうまく動かないかもしれません。

試行錯誤の履歴

まだまだ Pony を理解できていないので試行錯誤しながら作っています。

すんなりうまくいきました、ということだけ残して失敗は書かない方が書き手としては楽なのですが、読み手としては失敗を知る方が勉強になることが多いものなので、敢えて失敗を記録しておきます。

mpf_t, mpz_t のラッパークラスのあるべき Reference Capability がわからなかった

上にも書いたとおり、当初はデフォルトの ref にしていたのですが、ref だとアクターをまたがった受け渡しができません。

ref → val / val → ref の変換を行うようなコンストラクタなりメソッドなりを用意すれば逃げられるだろうとも思ったのですが、結局 val に変更しました。

これに伴い、mpz_t のビット操作関数はインプレースでビット操作を行う形から、ビット操作結果を返す形に変更しました（val の中身が変わるのはおかしいので）。

gmp_snprintf の結果をどうやって受ければいいかわからなかった

最初に以下のように書いてみたのですが、これではうまく動作しません：

// gmp_snprintf の結果を受け取る領域を Pony の文字列 s として作っておく
let s: String iso = recover iso String(bufSize) end
// gmp_snprintf の結果を Pony の文字列領域に書いてもらう
@__gmp_snprintf(s.cstring(), bufSize, "%Ff", _f)
// s に gmp_snprintf の結果が入っているはず？

cstring() は FFI 用にコピーした領域を返しているからです。String が保持している領域のポインタを返すだけの cpointer() を使えばよいのかとも思ったのですが、それでもだめでした。

そこで、malloc/free を使って以下のようにしてみました：

// gmp_snprintf の結果を受け取る領域を malloc しておく
let p: Pointer[U8] = @malloc(bufSize)
if p.is_null() then
  error
end
// gmp_snprintf の結果を malloc で確保した領域に書いてもらう
@__gmp_snprintf(p, bufSize, "%Ff", _f)
// さらに Pony の文字列に変換
let s: String ref = String.from_cpointer(p, bufSize)
// malloc で確保した領域を解放
@free(p)
// s に gmp_snprintf の結果が入っているはず？

これでもまだ正しくありません。@malloc で確保した領域は Pony ランタイムの管理下にないうえに @free で即座に解放してしまっているためです。String.from_cpointer ではなく String.copy_cpointer で値のコピーが必要です：

// gmp_snprintf の結果を受け取る領域を malloc しておく
let p: Pointer[U8] = @malloc(bufSize)
if p.is_null() then
  error
end
// gmp_snprintf の結果を malloc で確保した領域に書いてもらう
@__gmp_snprintf(p, bufSize, "%Ff", _f)
// さらに Pony の文字列に変換
let s: String ref = String.copy_cpointer(p, bufSize)
// malloc で確保した領域を解放
@free(p)
// s に gmp_snprintf の結果が入っている

実を言うと、これにも満足していません。本当に欲しいのは String ref でなく String val です。

仕方がないので String iso にちまちま一文字ずつコピーして iso を val にキャストする、などということをしていますがもっとスマートな方法がありそうな気がします。

まとめ

GMP を用いた多倍長演算を可能にするライブラリ Pony 向けに作成しました。
作ったライブラリを Corral で利用可能にしてみました（あまり自信はないので「してみました」レベル）。
多倍長演算できるようになったので、Pony で円周率の計算をしてみました。