Fortran で WebAssembly 一部実数関数を使う

# WebAssembly から javascript 関数利用

flang から WebAssembly を出力することが出来るのですが、ランタイムが無いため数学関数やメモリの allocate、I/O 等が出来ません。
[Fortran で WebAssembly 令和元年版](https://qiita.com/cure_honey/items/dbc06a088d2a7f86fbbd)


しかし、一部の数学関数は javascript の Math クラスにある関数を割り当てることで利用できることが分かりました。以下の記事が参考になりました。

[Fortran実装のゼータ関数をWASMにして動かしたかった件](https://abrakatabura.hatenablog.com/entry/2018/06/13/204136)

(私の以前の糞記事 [Fortran で WebAssembly　～または地獄めぐり～](https://qiita.com/cure_honey/items/ecf91523b9a5c6dfa25c)が恥ずかしながら多少なりとも参考になったようです。今度は私がお世話になりました、有難うございます。）

公式文書：
[WebAssembly JavaScript API を使用する](https://developer.mozilla.org/ja/docs/WebAssembly/Using_the_JavaScript_API)

## 要点
wasm を fetch してインスタンス化するときに、javascript 関数を付け加えられるので wasm-ld リンカで、「足りない yo！」と叱られた関数を補えばよいです。足りないものが分かったら **--allow-undefined** オプションで強行突破します。

```javascript
･･･　
      fetch('./logsin.wasm').then(response => response.arrayBuffer()
      ).then(bytes => {
        return WebAssembly.instantiate(bytes,
                    {env: {
                           __fd_sin_1 : (x) => Math.sin(x),
                           __fd_cos_1 : (x) => Math.cos(x),
                           __fd_tan_1 : (x) => Math.tan(x),
                           __fd_asin_1 : (x) => Math.asin(x),
                           __fd_acos_1 : (x) => Math.acos(x),
                           __fd_atan_1 : (x) => Math.atan(x),
                           __fd_atan2_1 : (y, x) => Math.atan2(y, x),
                           __fd_sinh_1 : (x) => Math.sinh(x),
                           __fd_cosh_1 : (x) => Math.cosh(x),
                           __fd_tanh_1 : (x) => Math.tanh(x),
                           __fd_exp_1 : (x) => Math.exp(x),
                           __fd_log_1 : (x) => Math.log(x),
                           __fd_log10_1 : (x) => Math.log10(x),
                           __mth_i_idnint : (x) => Math.round(x),
                           __mth_i_dhypot : (x, y) => Math.hypot(x, y),
                           __mth_i_dasinh : (x) => Math.asinh(x),
                           __mth_i_dacosh : (x) => Math.acosh(x),
                           __mth_i_datanh : (x) => Math.atanh(x)
                    }     }
            )
      }).then( r => { 云々
```
## 問題点
javascript は倍精度実数しかないので、関数も倍精度関数しかありません。特に問題となるのは複素関数で、上記参考記事にもありますが、構造体として戻り値をヒープに積むので、javascript 側からアクセスが難しいです。色々試して見ましたが、結局解決法が見つかりませんでした。

また複素関数のみならず、最適化をかけた場合に同じ引数に対する sin(x)、cos(x) が現れた場合、sin(x)、cos(x) を同時に計算して構造体として二値を返すルーチン (__fd_sincos_1(address, x)) が呼ばれるので、上記要点の素朴な方法では解決できません。今の所、最適化抑止させるしかありません。

（数値計算的には、同じ引数に対する sin(x)、cos(x) は、個々に計算するより、一緒に計算する方が必要な計算量が少なくなります。）

以下の記事によると、ECMAScript Module Integration Proposal でこのような問題に対する解決手段が与えられるようなので、無理せず機が熟するのを待てばいいのかもしれません。
[2019年のWebAssembly事情](https://qiita.com/bellbind/items/2619f8b71c3a69cc28be)

# 実例
## その1 最適化水準の問題

ここで　sincos = sin(x) * cos(x) と書いてしまうと、最適化水準に依らず __fd_sincos_1 関数が呼ばれてしまうので、冗長な表現で書いています。

なお javascrupt から呼び出す関数は、引数を value 属性にし、bind で c 互換名をつけておく必要があります。
### ソースプログラム
```fortran
module m_test
    implicit none
    integer, parameter :: kd = kind(0.0d0)
    real(kd), parameter :: pi = 4 * atan(1.0_kd)
contains

    real(kd) function sincos(x) bind(c, name = 'sincos')
        real(kd), value :: x
!        sincos = sin(x) * cos(x)
        real(kd) :: tmp
        tmp = sin(x)
        sincos = tmp * cos(x)
    end function sincos

    real(kd) function sin2(x) bind(c, name = 'sin2')
        real(kd), value :: x
        sin2 = sin(2 * x)
    end function sin2

end module m_test
```
###　コンソール入力
flang の最適化オプション -Oz をつけないことで、最適化を抑止し __fd_sincos_1 関数を呼ばずに、個々に __fd_sin_1, __fd_cos_1 関数を呼ばせます。

リンカ wasd-ld で数値組み込み関数が足りないと叱られるので、--allow-undefined オプションで強行突破します。足りない分は javascript で書いて、上記の要点に書いたように wasm のインスタンス化の時に与えます。

#### 最適化オプションをつけない場合

```
pc@VM10:~/WebAssembly$ flang --target=wasm32 -emit-llvm -c -S sincos.f90
warning: overriding the module target triple with wasm32 [-Woverride-module]
1 warning generated.
pc@VM10:~/WebAssembly$ llc --march=wasm32 -filetype=obj sincos.ll
pc@VM10:~/WebAssembly$ wasm-ld --no-entry --export-all  --strip-all -o sincos.wasm sincos.o
wasm-ld: error: sincos.o: undefined symbol: __fd_sin_1
wasm-ld: error: sincos.o: undefined symbol: __fd_cos_1
pc@VM10:~/WebAssembly$ wasm-ld --no-entry --export-all --allow-undefined --strip-all -o sincos.wasm sincos.o
pc@VM10:~/WebAssembly$
```

#### 最適化オプション -Oz をつけた場合
構造体をヒープに積んで返す関数 __fd_sincos_1 が要求されてしまいます。

```
pc@VM10:~/WebAssembly$ flang --target=wasm32 -emit-llvm -c -S -Oz sincos.f90
warning: overriding the module target triple with wasm32 [-Woverride-module]
1 warning generated.
pc@VM10:~/WebAssembly$ llc --march=wasm32 -filetype=obj sincos.ll
pc@VM10:~/WebAssembly$ wasm-ld --no-entry --export-all  --strip-all -o sincos.wasm sincos.o
wasm-ld: error: sincos.o: undefined symbol: __fd_sincos_1
wasm-ld: error: sincos.o: undefined symbol: __fd_sin_1
pc@VM10:~/WebAssembly$
```

### HTML file
```html
<!-- public/index.html -->
<!DOCTYPE html>
<html lang="en">
  <body>
    <script>

      fetch('./sincos.wasm').then(response => response.arrayBuffer()
      ).then(bytes => {
        return WebAssembly.instantiate(bytes,

                    {env: {
                           __fd_sin_1 : (x) => Math.sin(x),
                           __fd_cos_1 : (x) => Math.cos(x),
                           __fd_tan_1 : (x) => Math.tan(x),
                           __fd_asin_1 : (x) => Math.asin(x),
                           __fd_acos_1 : (x) => Math.acos(x),
                           __fd_atan_1 : (x) => Math.atan(x),
                           __fd_atan2_1 : (y, x) => Math.atan2(y, x),
                           __fd_sinh_1 : (x) => Math.sinh(x),
                           __fd_cosh_1 : (x) => Math.cosh(x),
                           __fd_tanh_1 : (x) => Math.tanh(x),
                           __fd_exp_1 : (x) => Math.exp(x),
                           __fd_log_1 : (x) => Math.log(x),
                           __fd_log10_1 : (x) => Math.log10(x),
                           __mth_i_idnint : (x) => Math.round(x),
                           __mth_i_dhypot : (x, y) => Math.hypot(x, y),
                           __mth_i_dasinh : (x) => Math.asinh(x),
                           __mth_i_dacosh : (x) => Math.acosh(x),
                           __mth_i_datanh : (x) => Math.atanh(x)
                    }     }
            )
      }).then( r => {
        const PI = Math.PI
        const ex = r.instance.exports
        const res1 = ex.sincos(PI / 6) * 2
        const res2 = ex.sin2(PI / 6)
        document.querySelector("#result1").innerHTML = res1
        document.querySelector("#result2").innerHTML = res2
      })
    </script>

    sin(x)cos(x)
    <div id="result1"></div>
    sin(2x)
    <div id="result2"></div>

  </body>
</html>
```

```
pc@VM10:~/WebAssembly$ python -m SimpleHTTPServer 8080
Serving HTTP on 0.0.0.0 port 8080 ...
```

#### 出力結果
sin(2x) = 2 * sin(x) * cos(x) を確認しています。(x=pi/6)

```
sin(x)cos(x)
0.8660254037844386
sin(2x)
0.8660254037844386
```
![sincos.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44522/4fb81ddc-239d-8b06-94e8-cab7293fdc6f.png)

## その2 log(sin(x)) = -log(2) - $\sum_{i=1}^\infty{\cos(2ix)\over i}$ 

### ソースプログラム
``` fortran
module logsin_m
    implicit none
    integer, parameter :: kd = kind(0.0d0)

contains

    real(kd) function cossum(n, x) result(s)
         integer , intent(in) :: n
         real(kd), intent(in) :: x
         integer :: i
         s = 0.0_kd
         do i = 1, n
             s = s + cos(2 * i * x) / i
         end do
    end function cossum

    real(kd) function alogsin1(n, x) result(res) bind(c, name = 'alogsin1')
        integer , value :: n
        real(kd), value :: x
        res = -log(2.0_kd) - cossum(n, x)
    end function alogsin1

    real(kd) function alogsin2(x) result(res) bind(c, name = 'alogsin2')
        real(kd), value :: x
        res = log(sin(x))
    end function alogsin2

end module logsin_m
```

### コンソール入力
```
pc@VM10:~/WebAssembly$ flang --target=wasm32 -emit-llvm -c -S -Oz logsin.f90
warning: overriding the module target triple with wasm32 [-Woverride-module]
1 warning generated.
pc@VM10:~/WebAssembly$ llc --march=wasm32 -filetype=obj logsin.ll
pc@VM10:~/WebAssembly$ wasm-ld --no-entry --export-all --allow-undefined --strip-all -o logsin.wasm logsin.o
```

### HTML file
```html
<!-- public/index.html -->
<!DOCTYPE html>
<html lang="en">
  <body>
    <script>

      fetch('./logsin.wasm').then(response => response.arrayBuffer()
      ).then(bytes => {
        return WebAssembly.instantiate(bytes,

                    {env: {
                           __fd_sin_1 : (x) => Math.sin(x),
                           __fd_cos_1 : (x) => Math.cos(x),
                           __fd_tan_1 : (x) => Math.tan(x),
                           __fd_asin_1 : (x) => Math.asin(x),
                           __fd_acos_1 : (x) => Math.acos(x),
                           __fd_atan_1 : (x) => Math.atan(x),
                           __fd_atan2_1 : (y, x) => Math.atan2(y, x),
                           __fd_sinh_1 : (x) => Math.sinh(x),
                           __fd_cosh_1 : (x) => Math.cosh(x),
                           __fd_tanh_1 : (x) => Math.tanh(x),
                           __fd_exp_1 : (x) => Math.exp(x),
                           __fd_log_1 : (x) => Math.log(x),
                           __fd_log10_1 : (x) => Math.log10(x),
                           __mth_i_idnint : (x) => Math.round(x),
                           __mth_i_dhypot : (x, y) => Math.hypot(x, y),
                           __mth_i_dasinh : (x) => Math.asinh(x),
                           __mth_i_dacosh : (x) => Math.acosh(x),
                           __mth_i_datanh : (x) => Math.atanh(x)
                    }     }
            )
      }).then( r => {
        const PI = Math.PI
        const ex = r.instance.exports
        const x = PI/6
        const res1 = ex.alogsin1(10e5, x)
        const res2 = ex.alogsin2(x)
        document.querySelector("#x").innerHTML = x
        document.querySelector("#result1").innerHTML = res1
        document.querySelector("#result2").innerHTML = res2
      })
    </script>
    x
    <div id="x"></div>
    -log(2)-\sum_i^n \cos(2 i x) / i
    <div id="result1"></div>
    log(sin(x))
    <div id="result2"></div>

  </body>
</html>
```

### 出力結果
和を十万項 (10^5) まで取った場合の比較。10^6 まで計算すると瞬時には結果が出ず一瞬かかります。式から見て収束は遅そうで、小数点以下 5 桁目あたりまでは一致しそうです。
$$\log(\sin(x)) = -\log(2) - \sum_{i=1}^\infty{\cos(2ix)\over i}$$ 

x=pi/6 を与えました。sin(pi/6)=0.5 なので log(0.5)=-0.69314718.. が理論値。

```
x
0.5235987755982988
-log(2)-\sum_i^n \cos(2 i x) / i
-0.6931461805604385
log(sin(x))
-0.6931471805599454
```
果たしておおむね予想通りの結果の精度まで一致しています。

![logsin.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44522/1d262760-8248-a7fd-7d1c-091405fd69a0.png)