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

背景

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

#julialang でGR使ってゼータ関数をプロット。クリップの仕方が分からなかったから関数を再定義して調整。 https://t.co/qdWewxlXow pic.twitter.com/U8vXMEGzDl
— kjunichi (@kjunichi) 2018年5月14日

#Julia言語

私もやってみました。https://t.co/g94D7R57BT

zetafは

zetaf(s) = min(abs(zeta(s)), 3)

の1行で定義可能だし, zは

z = zetaf.(x' .+ im*y)

で作れます。Julia言語では . をうまく使うとコードがシンプルになります。続くhttps://t.co/CNglo4dOX3
— 黒木玄 Gen Kuroki (@genkuroki) 2018年5月14日

そんな訳でゼータ関数をもっと身近に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でなくてもWindowsかLinuxなら普通にインストールできると思われる。

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

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

こうして、実数の定義域の三角関数や指数関数の外部参照のみとなったので、あれ？直接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でビルドして動かして、計算結果のプロットを試してみた。

他のサイトで表示されるようなグラフなのでまぁ、だいたい良いのではなかろうか。

おまけその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))

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

その２に続く。。。

ゼータ関数つながり

14年前の記事

Nikkei Linux 2004年7月号

4年前の記事

Castoを知りkqueueにたどり着いた件

non vorrei lavorare

昔はおもにプログラミングやガジェット系、今は？

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

背景