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

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

# はじめに

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

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

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

# 方針

[チュートリアル](https://tutorial.ponylang.io/expressions/ops.html) にもあるように、Pony では `+`, `-`, `*`, `/`,... のオペレータをメソッドとして定義できます。

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

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

```pony
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 の結果を返す
  (略)
```

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

```pony
  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)
```

このコードは [円周率を求めるプログラムを書く](https://qiita.com/PenguinCabinet/items/d9ed9aacec4149e72efc) の記事を見て書きました。見比べて頂ければそのままなのがわかると思います。

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

```pony
  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 のドキュメントにある

```c++
Function: void mpf_init (mpf_t x)
```

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

```c++
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 用パッケージ作成

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

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

1. [Library Scaffolding Generator](https://github.com/ponylang/library-scaffolding-generator/) を使ってパッケージの雛形を作成
2. 作った雛形に、作成した GMP ラッパーソースを配置
3. corral.json にパッケージ情報を記述すると共に、post_fetch_or_update で make all を実行するように記述（共有ライブラリを自動的にビルドするため）

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

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

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

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

# 試行錯誤の履歴

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

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

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

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

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

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

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

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

```pony
// 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` を使って以下のようにしてみました：

```pony
// 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` で値のコピーが必要です：

```pony
// 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 を用いた多倍長演算を可能にするライブラリ](https://github.com/tadashi9e/gmp4pony) Pony 向けに作成しました。
* 作ったライブラリを Corral で利用可能にしてみました（あまり自信はないので「してみました」レベル）。
* 多倍長演算できるようになったので、[Pony で円周率の計算](https://github.com/tadashi9e/gmp4pony/blob/main/examples/pi/pi.pony) をしてみました。