non vorrei lavorare

ブログ名の通りです。javascript three.js mruby rust OCaml golang julialang blender

Fortran実装のゼータ関数をWASMにして動かしたかった件

こんばんは、5歳の誕生日を過ぎ、ついに次男がよるもパンツで寝るようになりました!おむつ代が家計を苦しめていたので助かりますw

背景

Juliaでゼータ関数をプロットしたインスタをツイートしたら本職の方に良いこと教えてもらったことが 最近あり、やっぱり何事も情報発信は大事だと強く思った出来事があった。

そんな訳でゼータ関数をもっと身近にWebブラウザ上に表示したいと考えた。

Fortran実装のゼータ関数のコードを見つけた

tsujimotter.hatenablog.com

Rubyの実装を知り、さらに

slpr.sakura.ne.jp

を見つけた。ここにFortranで実装したゼータ関数のコードがあった。 Ruby実装よりアルゴリズムレベルで高速な方式が採用されているとのことだったので、 これをはじめはCに移植して、Emscriptenコースを考えた。

Fortranをf2cで変換すればEmscriptenでWASMでは?

イマドキはf90らしい。f2cはf77が対象だった。

f90をf77に変換してcにする手も探したが。。。

NVIDIAがflangを作っていた

幸いgfortranでもf90はビルド出来た。 しかし、試しに付属のf2cを試したが、NG。

llvmベースのFortranのflangを知った。

github.com

Qiitaでズバリな神記事!が!!!

qiita.com

flangのインストール

www.scivision.co

ここで紹介されているようにcondaを使うと超簡単!

今回はmacOS上で作業を始めてしまったので、dockerを使ってUnbutu18のイメージを用意して、これにcondaを入れて、 そこにflangを入れた。

apt-get update
apt install wget
wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh 
source .bashrc
conda install -c conda-forge flang

こんな感じでインストールした模様。(historyコマンドから拾い出した)

別に、dockerでなくてもWindowsLinuxなら普通にインストールできると思われる。

Fortranのランタイム問題

flang経由で生成したWASTにはいくつかのFortranのランタイム関数を外部参照していた。

wastを生成すると、以下のような外部参照シンボルが明らかに。

 (import "env" "__mth_i_cdabs" (func $__mth_i_cdabs (param f64 f64) (result f64)))
 (import "env" "__mth_i_cddiv" (func $__mth_i_cddiv (param i32 f64 f64 f64 f64 i32) (result i32)))
 (import "env" "__mth_i_cdexp" (func $__mth_i_cdexp (param i32 f64 f64 i32) (result i32)))
 (import "env" "__mth_i_cdlog" (func $__mth_i_cdlog (param i32 f64 f64 i32) (result i32)))
 (import "env" "__mth_i_cdpowcd" (func $__mth_i_cdpowcd (param i32 f64 f64 f64 f64 i32) (result i32)))
 (import "env" "__mth_i_cdsin" (func $__mth_i_cdsin (param i32 f64 f64 i32) (result i32)))
 (import "env" "__mth_i_dlog" (func $__mth_i_dlog (param f64) (result f64)))
 (import "env" "__mth_i_idnint" (func $__mth_i_idnint (param f64 i32) (result i32)))
(import "env" "f90_mzeroz16" (func $f90_mzeroz16 (param i32 i64 i32)))

これは、flangのコードを調べたら、C実装であった。 そこで、Emscriptenでこの実装を用意して、WASTを組み合わせれば、参照問題を解決できるのでは? と考え、これらの関数を片っ端からCで実装した。

GitHubのflangのリポジトリのコードを手動でマクロ展開して コードに落とすだけだが。一部、アセンブリで高速化しているせいか、引数がCの定義とWASTで異なっているものも あったが、WASTに寄せて対応した。)

モリーエラー発生

こうしてFortranとCを合わせたWASTを何とかWASMに変換して 動かしたが、メモリーエラー。

さらりとFortranとCのwastを合せたと書いたが、flangのllvmが結構古く、これを元にしたwastの形式も 微妙にEmscriptenのそれと異なり、この辺りを手動で吸収して対応した。

モリーエラー原因は、関数呼び出し時に構造体があると、WASMはこれをグローバルアドレスで処理しようとして、 Fortran側のコードのアドレスと、C側のコードのアドレスが不一致であった為の模様。

実際、シンプルなコードで一致させたら呼び出せた。

しかし、コンパイラで吐かれたFortranとCのコードで全ての箇所を対応するには気が遠くなった。

実はEmscripten要らなかった

構造体を使うのは困難な事が分かったので、これを解決すべく、Fortran側で、構造体を使わないコードに 展開するという作業を行った。

今回の場合で言えば、Fortranで用意されている複素数関数をすべて、実数ベースの関数に展開するということに 相当した。

今回ラップした関数たちは以下の通り。

function mylog(s)
  implicit none
  complex(kind(0d0)),intent(in)::s
  complex(kind(0d0))::mylog
  double precision x,y

  x = atan2(aimag(s),real(s))
  y = log(hypot(real(s), aimag(s)))

  mylog = cmplx(x,y)
end function

function mysin(s)
  implicit none
  complex(kind(0d0)),intent(in)::s
  complex(kind(0d0))::mysin
  double precision x1,y1,x,y

  x1 = real(s)
  y1 = aimag(s)

  x = sin(x1)
  y = cos(x1)
  x = x * cosh(y1)
  y = y * sinh(y1)

  mysin = cmplx(x,y)
end function

function mydiv(s1,s2)
  implicit none
  complex(kind(0d0)),intent(in)::s1,s2
  complex(kind(0d0))::mydiv
  double precision::x1,y1,x2,y2,x,y,r,d,r_mag,i_mag


  x1 = real(s1)
  y1 = aimag(s1)
  x2 = real(s2)
  y2 = aimag(s2)

  r_mag = x2
  if(r_mag <0) then
          r_mag = -r_mag
  end if
  i_mag = y2
  if(i_mag < 0) then
          i_mag = -i_mag
  end if
  if(r_mag <= i_mag) then
          r = x2/y2
          d = 1.0 /(y2 * (1 + r*r))
          x = (x1 * r + y1) * d
          y = (y1 * r + x1) * d
  else
          r = y2/x2
          d = 1.0/(x2 * (1 + r*r))
          x = (x1 + y1 * r) * d
          y = (y1 - x1 * r) * d
  end if
  mydiv = cmplx(x,y)
end function

function mydivr1(r1, s2)
  implicit none
  double precision,intent(in)::r1
  complex(kind(0d0)),intent(in)::s2
  complex(kind(0d0))::mydivr1
  complex(kind(0d0))::s1

  s1 = cmplx(r1,0.0)
  mydivr1 = mydiv(s1,s2)
end function

function mydivr2(s1, r2)
  implicit none
  double precision,intent(in)::r2
  complex(kind(0d0)),intent(in)::s1
  complex(kind(0d0))::mydivr2
  complex(kind(0d0))::s2

  s2 = cmplx(r2,0.0)
  mydivr2 = mydiv(s1,s2)
end function

function mypow(s1,s2)
  implicit none
  complex(kind(0d0)),intent(in)::s1,s2
  complex(kind(0d0))::mypow
  double precision::x1,y1,x2,y2,logr,logi,x,y,z,w

  x1 = real(s1)
  y1 = aimag(s1)
  x2 = real(s2)
  y2 = aimag(s2)

  logr = log(hypot(x1,y1))
  logi = atan2(y1,x1)

  x = exp(logr * x2 - logi * y2)
  y = logr * y2 + logi * x2

  z = x * cos(y)
  w = x * sin(y)

  mypow = cmplx(z,w)  
end function mypow

function mypowr1(r1,s2)
  implicit none
  double precision,intent(in)::r1
  complex(kind(0d0)),intent(in)::s2
  complex(kind(0d0))::mypowr1
  complex(kind(0d0))::s1

  s1 = cmplx(r1,0)
  mypowr1 = mypow(s1,s2)
end function mypowr1

function myexp(s)
  implicit none
  complex(kind(0d0)),intent(in)::s
  complex(kind(0d0))::myexp
  double precision::x,y,a,b

  x = real(s)
  y = aimag(s)

  a = exp(x)*cos(y)
  b = exp(x)*sin(y)

  myexp = cmplx(a,b)
  
end function myexp

このおかげで、複素数関数に強くなったかもw。

こうして、実数の定義域の三角関数や指数関数の外部参照のみとなったので、あれ? 直接JavaScriptのMath.xxxを叩けば良いんじゃね?

複素数関数を展開後のWASTでの外部参照の様子

 (import "env" "__mth_i_cdabs" (func $__mth_i_cdabs (param f64 f64) (result f64)))
 (import "env" "__mth_i_datan2" (func $__mth_i_datan2 (param f64 f64) (result f64)))
 (import "env" "__mth_i_dcos" (func $__mth_i_dcos (param f64) (result f64)))
 (import "env" "__mth_i_dcosh" (func $__mth_i_dcosh (param f64) (result f64)))
 (import "env" "__mth_i_dexp" (func $__mth_i_dexp (param f64) (result f64)))
 (import "env" "__mth_i_dhypot" (func $__mth_i_dhypot (param f64 f64 i32) (result f64)))
 (import "env" "__mth_i_dlog" (func $__mth_i_dlog (param f64) (result f64)))
 (import "env" "__mth_i_dsin" (func $__mth_i_dsin (param f64) (result f64)))
 (import "env" "__mth_i_dsinh" (func $__mth_i_dsinh (param f64) (result f64)))
 (import "env" "__mth_i_idnint" (func $__mth_i_idnint (param f64 i32) (result i32)))

こうして、試行錯誤して出来たEmscriptenフリーなWASMを動かすためのJavaScriptのコード がこちら (Mozillaが最近公開したWASMのIDE環境のmain.jsを大いに参考にしてつくられている)

//const fetch = require('node-fetch');
const fs = require('fs')
let x = './zeta4.wasm'

let instance = null

let s = "";
fs.readFile(x, (err,bytes)=>{
  WebAssembly.instantiate(bytes, {
    env: {
      __mth_i_idnint: function(d) {
        console.log(`d = ${d}`)
        if(d>0) {
          return (d+0.5)|0
        } else {
          return (d-0.5)|0
        }
      },
      __mth_i_dsin: function(d) {
        console.log(`d = ${d}`)
        return Math.sin(d)
      },
      __mth_i_dsinh: function (d) {
        console.log(`d = ${d}`)
        return Math.sinh(d)
      },
      __mth_i_dlog: function(d) {
        return Math.log(d)
      },
      __mth_i_dexp: function (d) {
        console.log(`d = ${d}`)
        return Math.exp(d)
      },
      __mth_i_dhypot: function(x,y,d) {
        return Math.hypot(x,y)
      },
       __mth_i_dcos: function (d) {
        console.log(`d = ${d}`)
        return Math.cos(d)
      },
      __mth_i_dcosh: function (d) {
        console.log(`d = ${d}`)
        return Math.cosh(d)
      },
      __mth_i_datan2: function (x,y) {
        console.log(`x,y = ${x},${y}`)
        return Math.atan2(x,y)
      },
      __mth_i_cdabs: function (x,y) {
        return Math.sqrt(x*x + y*y)
      },
      putc_js: function (c) {
        c = String.fromCharCode(c);
        if (c == "\n") {
          console.log(s);
          s = "";
        } else {
          s += c;
        }
      }
    }
  }).then(results => {
    //console.log(results)
    instance = results.instance
    console.log( instance.exports.zeta3(0.5,14.0))
  }).catch(console.error) 
})

再度メモリーエラー発生

こうして、紆余曲折して、flangから生成したllファイルのみから生成したWASMを動かした!

しかし!

モリーエラー!

メモリ操作が鬼門、複素三角関数や指数関数はクリアできたが、 神Qiita記事にもあるようにFortranの配列処理が駄目らしい 記事によると、サイズ固定ならいけそうにも解釈できるのだけど。。

学べたこと

関数の引数が構造体だとFortranでWASTにしたものと、EmscriptenでWASTにしたものとで、 そのままでは、アドレス空間が違うので、WASM実行時にエラーになる。 (WAST->WASMは変換出来てしまう)

前述の問題があり、これを簡単に避ける方法としては、 定義域が複素数なFotrtanな関数はFortran側で定義域が実数の関数での計算に展開しておく 必要がある。

wasmをCに変換するツールを試したが、機械的に翻訳するだけのようで、wastの1文を一対一にそのまま Cに変換するだけだった。そのため、構造体の受け渡し等もWASTでのグローバルアドレス値を渡すということいに なり、とてもFortranをCに変換するといった用途には使えないことが分かった。

おまけ

今回対象にしたFortran実装のゼータ関数はflangでの実績は明記されておらず、 flangで動くのかを検証すべく、 フツーにflangでビルドして動かして、計算結果のプロットを試してみた。

-1 0 1 2 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 y1 他のサイトで表示されるようなグラフなのでまぁ、だいたい良いのではなかろうか。

おまけその2

今回の作業の途中でJavaScript実装のゼータ関数を見つけた。 こちらのSinc関数の回転体と比べると分かるのだが、かなり 遅い、というか、滅茶滅茶遅いので、WebWorkerを駆使してマルチコア対応もしてが、 それでも遅い。。

See the Pen abs(Zeta) by Junichi Kajiwara (@kjunichi) on CodePen.

なんとか早さを売りにしている今回見つけた Fortan実装のゼータ関数を動かせたら、さくっと表示されるのかなぁと夢見て 作業をした。

Sinc関数の回転体

See the Pen rot sinc by Junichi Kajiwara (@kjunichi) on CodePen.

PS

どうやら

RuntimeError: memory access out of bounds

は配列操作が原因ではなかった。。。

wastに

(return (f64.const 1))

を頭の方から書いてはずらして、落ちる箇所を突き止めたのだった。

その2に続く。。。

関連記事

14年前の記事

4年前の記事