自作プログラミング言語Lazeで超級積分を解いてみた話

まえがき

（友達が）自作プログラミング言語Lazeを作りました。その言語で積分ガチャの超級積分を解いてみようという話です。ちなみにオープンソースで開発されてて、公式ウェブサイトもあります。興味があればどうぞ（ダイマ）

積分を解こう

積分を解くといってもどこぞの計算知能みたいに問題を入れれば解いてくれるみたいなマシーンではありません。タイトル詐欺でゴメンナサイ。区分求積法という方法を用いて定積分の値を小数で求めようという算段です。積分力に自信のある人は一緒に解きながら読むと楽しいかもしれません。

プログラム

// 積分する関数
関数:f (実数:x) => (実数:y) {
    y = x * x;
}

// 区分求積
関数:定積分 (実数: 下端, 実数: 上端) => (実数: 結果) {
    整数:n = 10000;
    結果 = 0.0;
    (整数:k = 下端 * n;) から (k >= (上端 * n)) まで (k += 1;) {
        結果 += f(1.0 * k/n);
    }
    結果 = 結果 / n;
}

// エントリポイント
関数:実行 () => () {
    実数:下端 = 0.0;
    実数:上端 = 1.0;
    表示(定積分(下端, 上端));
}

解説

Lazeは母国語でプログラミングできる言語なので、日本人でも結構読みやすい気がしますが一応解説。

2～4行目で被積分関数f(x)を定義しています。
7～14行目で区分求積法を実行してます。1/10000ずつ関数の値を足し算していくことで関数の面積を近似しています。区分求積法は関数を細かく切って足し合わせて面積を近似する(n→∞で一致する)ことです。

参考＆図の出典：おいしい数学

最後に17～21行目で範囲を決めて定積分を行っています。

いざ積分

長々と解説してしまいましたがいざプログラムを動かしてみましょう。

第一問

答え

数学がメインではないので解法はばっさりカットしました。これをLazeで求めてみましょう。

---------- START ----------
0.287079478363991

答えは少数に直すと0.2870465829187831となるのでかなりの精度(誤差0.011%)で近似できていることがわかります。n=10000でこの値というだけでn=100000では0.2870498722725148(誤差0.0011%)、n=1000000では0.28704691185224707(誤差0.00011%)とnの値を大きくしていくと段々近似されていきます。

第二問

答え

不定積分が現れました。このままだと解けない（区分求積法で値が出せない）ので適当に範囲を設定してやる必要があります。範囲が大きすぎると時間がかかるので今回は適当に0～1にしました。Lazeを使って値を出してみると...

---------- START ----------
0.11438084084341338

(n=10000のとき)

f(1) - f(0) = 0.1143957049107332なのでしっかりと近似できています。

第三問

答え

ものすごい複雑な形ですがLazeに計算できるのでしょうか？範囲は1～2(0だと発散する)で計算してみます。

---------- START ----------
0.5826849005818652

f(2) - f(1) = 0.5822617392595333となって十分に値が求められていることが分かります。

あとがき

自作プログラミング言語Lazeを使って超級積分たちを解いてみました。結構プログラミング要素が少なくなった感もありますが初投稿なので大目に見て下さい💦Qiitaを読んでいるような人たちはプログラミングができる人が多いと思いますが、周りにプログラミングを始めようと思ってるけど英語が難しいな～とか環境構築ができない～みたいな人に勧めて見て下さい。私が喜びます。

おまけ

プログラミング言語には通常あたりまえのように三角関数やら指数関数やら対数関数やら平方根やらが標準ライブラリ(標準機能)で用意されています。これらの内部実装を気にしたことがある人は少ないと思いますが、Lazeは自作言語なのでもちろん自分の手で実装する必要があります。使えるものは四則演算のみです。さてどうやってやるんだというのがおまけです。

三角関数

sin, cos, tanの3つがありますがtanはsinとcosがわかれば求められるので実質sinとcosの2つを求めればいいことになります。ここで出てくるのがマクローリン展開です。x=0近辺で多項式だけでいろんな関数を近似できるという優れものです。

出典：https://manabitimes.jp/math/570

x=0に近いほど精度が高いわけですが、三角関数の周期性とsinとcosの性質を利用することで0～π/4までの範囲だけの近似で済みます。

実際の実装

関数:sin (実数:input) => (実数:結果) {
  実数:PI = 3.14159265358979323846264338;
  整数:temp = input / (2.0 * PI);
  実数:rad = input - (2.0 * PI * temp);
  実数:plusminus = 1.0;
  もし (rad < 0.0) ならば {
    rad = rad + (2.0 * PI);
  }
  もし (rad > PI) ならば {
    plusminus = -1.0;
    もし (rad > (1.5 * PI)) ならば {
      rad = (2.0 * PI) - rad;
    } でなければ {
      rad = rad - PI;
    }
  } でなければ {
    もし (rad > (0.5 * PI)) ならば {
      rad = PI - rad;
    }
  }
  もし (rad > (0.25 * PI)) ならば {
    rad = (0.5 * PI) - rad;
    実数:doubleRad=rad*rad;
    結果 = 1.0 - (doubleRad)/2.0 + (doubleRad*doubleRad)/24.0 - (doubleRad*doubleRad*doubleRad)/720.0 + (doubleRad*doubleRad*doubleRad*doubleRad)/40320.0 - (doubleRad*doubleRad*doubleRad*doubleRad*doubleRad)/3628800.0 +(doubleRad*doubleRad*doubleRad*doubleRad*doubleRad*doubleRad)/479001600.0;
    結果 = 結果 * plusminus;
  } でなければ {
    実数:doubleRad=rad*rad;
    結果 = rad - (rad*doubleRad)/6.0 + (rad*doubleRad*doubleRad)/120.0 - (rad*doubleRad*doubleRad*doubleRad)/5040.0 + (rad*doubleRad*doubleRad*doubleRad*doubleRad)/362880.0 - (rad*doubleRad*doubleRad*doubleRad*doubleRad*doubleRad)/39916800.0 + (rad*doubleRad*doubleRad*doubleRad*doubleRad*doubleRad*doubleRad)/6227020800.0 - (rad*doubleRad*doubleRad*doubleRad*doubleRad*doubleRad*doubleRad*doubleRad)/1307674368000.0;
    結果= 結果 * plusminus;
  }
}

鬼長い部分がマクローリン展開を実装している部分ですね

指数関数

三角関数同様にマクローリン展開展開して求めます。ただしこれに関しては周期性はないので指数によってある程度高い次数まで展開する必要があるかもしれません。

出典：https://manabitimes.jp/math/570

実際の実装

関数: eˣ = (実数:x) => (実数:y) {
  y = 1.0 + x + ((1.0/2.0)*x*x) + ((1.0/6.0)*x*x*x) + ((1.0/24.0)*x*x*x*x) + ((1.0/120.0)*x*x*x*x*x) + ((1.0/720)*x*x*x*x*x*x);
}

今回は6乗までの近似にしましたがこの程度の近似では、ある程度正確な値が求まるのはx<1の範囲ぐらいです。今回の積分では0<x<1の範囲だったので大きな問題はなかったのですが、もっと大きいxを代入するときは別の近似方法が必要かもしれません。

対数関数

対数関数は-1<x<1の範囲でしかマクローリン展開が行えません。逆数をとることでマクローリン展開に持ち込むことも可能ですが、今回は別のアプローチで挑もうと思います。

1/xの積分がlog|x|になることを利用して1/xを0～aまでの範囲で区分求積することでlog|a|の値を求めることができます。nの値を大きくすれば簡単に精度を上げられるので比較的お手軽だと思います。ただaの値が大きくなってくると処理の時間の割に精度が落ちるのでその場合は逆数を取る方法を取った方がいいと思います。

やってることはただの区分求積なので実装は省きます。

平方根

平方根を求めるプログラムは2乗すれば元の値になることを利用すれば、頑張れば初心者でも書けなくはないのですが、方法は計算量が非常に多いことが問題になるので数学的な法則を利用して求めます。

この漸化式を利用することで√aを求めることができます。この漸化式は非常に収束が早く、初項の値にほとんど関係なく5回程度で概ね収束します。

実際の実装

関数:root (実数:x) => (実数:y) {
  y = 5;
  (整数:i = 0;) から (i == 20) まで (i++;) {
    y = 0.5 * (y + (x / y));
  }
}

この場合では初項5で20回演算を行っています。doubleの範囲(64bit)ではほぼ一致する程度の精度が出せているはずです。

参考

マクローリン展開：https://manabitimes.jp/math/570