数値計算する（シェルスクリプト）

シェルスクリプト

課題

なんだか無性に円周率を求めたくなった。数値計算で求めたい。

解決策

バーゼル問題の等式を利用する（以下）。

\displaystyle \sum_{n=1}^\infty \frac{1}{n^2} = \frac{\pi^2}{6}


数値計算をイメージしやすくするために、最初の一部の項をベタ書きで展開してみる。これを参考に数値計算のスクリプトを組み上げる。

\displaystyle \sum_{n=1}^6 \frac{1}{n^2} = \frac{1}{1^2} + \frac{1}{2^2} + \frac{1}{3^2} + \frac{1}{4^2} + \frac{1}{5^2} + \frac{1}{6^2} \approx \frac{\pi^2}{6}


①整数列を作る

$ seq 1 6
1
2
3
4
5
6

②二乗する

$ seq 1 6 | awk '{print $1*$1}'
1
4
9
16
25
36

③逆数を取る

$ seq 1 6 | awk '{print $1*$1}' | awk '{print 1/$1}'
1
0.25
0.111111
0.0625
0.04
0.0277778

④合計する

$ seq 1 6 | awk '{print $1*$1}' | awk '{print 1/$1}' | awk '{sum+=$1} END {print sum}'
1.49139

⑤6倍する

$ seq 1 6 | awk '{print $1*$1}' | awk '{print 1/$1}' | awk '{sum+=$1} END {print sum}' | awk '{print $1*6}'
8.94834

⑥平方根をとる

$ seq 1 6 | awk '{print $1*$1}' | awk '{print 1/$1}' | awk '{sum+=$1} END {print sum}' | awk '{print $1*6}' | awk '{print sqrt($1)}'
2.99138


これで円周率の値が数値計算できるはずだが、もとの等式は無限級数であるため、先頭6項のみで計算しても、円周率といえる値にはなっていない。この項数を大きくしたとき、どのような値になるか、いくつか試してみた。項数が増加するにつれ、円周率の真の値に近づいていることがわかる。

項数 値
6 2.99138
10 3.04936
100 3.13207
1000 3.14063
10000 3.14149

所感

この記事の目的はもちろん円周率を求めることではない（←！？）。単純な演算をパイプラインで複数接続することで複雑な演算を実現できるということを示すための良い題材だと思ったからだ。この例ではほとんどをawkで記述しているが、同じawkなら全部まとめればよいのでは？と思った人には、おそらくこの記事で言いたいことは伝わっていないだろう。