数学関数を入れ始める

VBAHaskellには数学関数がまだない。四則演算、指数関数、対数関数、多項式、min、max 程度が Haskell_2_stdFun にあるのと、乱数生成機が misc_random にあるだけだ。 とりあえず数学関数を少しずつ書いていくためのbasファイルを追加した。

VBA/misc_math.bas at master · mYmd/VBA · GitHub

中身はまだほとんどない。仕事で使うわけでもないからなかなか進捗しないだろう。
とりあえず Simpson法 による簡単な数値積分の関数を入れたので、そのことを書く。もちろん表面的にはループを使いたくない。

' シンプソン法による数値積分
Function integral_simpson(ByRef fun As Variant, _
                            ByVal begin_ As Double, _
                            ByVal end_ As Double, _
                            ByVal n As Long) As Double
    Dim xs As Variant, ys As Variant
    ' 独立変数の列
    xs = mapF(p_poly(, Array((end_ - begin_) / 2 / n, begin_)), iota(0, 2 * n))
    ' 関数の値を計算
    ys = mapF(fun, xs)
    ' 定数列 (1, 4, 2, 4, 2, 4, ... , 2, 4, 1) を作る
    Dim constants As Variant
    ReDim constants(0 To 2 * n)
    Call fillPattern(constants, Array(2, 4))
    constants(0) = 1
    constants(2 * n) = 1
    ' 積分値の計算
    integral_simpson = foldl1(p_plus, zipWith(p_mult, constants, ys)) * (end_ - begin_) / n / 6
End Function

4つの引数の意味は以下の通りで、返り値の型は Double だ。

引数	型	意味	備考
fun	ByRef Variant	対象関数	VBAHaskellで定義されている関数
begin_	ByVal Douuble	積分区間の下限
end_	ByVal Douuble	積分区間の上限
n	ByVal Long	分割数	それぞれの中点でも分割される

最後の備考にあるように、分割数を 100 とすると内部的には 200 分割されて、関数値は 201 個計算されることになる。Simpson法の仕様により、こうするのが自然に書ける。
念のため書いておくと、'constants' という配列は下図の「重み」に相当する部分だ。

円の面積を求めてチェックする。

 ' 扇型の面積を4倍する （ 区間 [0, 1] を1,000,000 に分割して求める ）
 Debug.Print 4 * integral_simpson(p_pow(p_minus(1, p_pow(,2)),0.5), 0,1,1000000)
 3.1415926534273 

 ' 真の値に近いもの
 Debug.Print 4*atn(1)
 3.14159265358979 

p_pow はべき乗を求める関数 pow_fun のファンクタで、p_pow(p_minus(1, p_pow(,2)),0.5) は √(1 - x ^ 2) を表すことになる。
100万分割で真の値と少数以下9桁まで一致した。1000分割でも4桁まで一致する。
表面的にはループはないもののスピードは遅い。100万分割だと3秒ほどかかってしまう。

中で1次元配列を他の1次元配列の繰り返しで埋める fillPattern という Sub を使っており、今回 Haskell_4_vector.bas に追加した。
ついでに多項式関数 poly もかなり遅かったので改善した。これは Haskell_2_stdFun にある。両方とも素直にループを使って実装している。このような汎用部品はパフォーマンスを優先するためしかたがない。