円周率を積分で近似計算する

はじめに

どうもこんにちは、🐟️🍊みかん🍊🐟️です。 こんなこと を言ってしまったので、今回は円周率などを高校範囲で有理数評価する方法を書いていきます。この方法は中学生の時に思い付いた方法ですが、いくつかの計算において一般的な計算方法がわかったために簡略化が加わったため、計算自体は非常に簡単なものになっています。そして、僕が知る限り高校範囲内において、一般的に多少の桁数を(高校範囲のみで)求める際には 計算が軽い ことと、他の数学定数の計算にも応用できる点でかなり使い勝手のいい方法だと思います。
本稿を書くにあたって、先行研究を全く調べていないので恐らく既知だと思います。また、$\pi$の$50$桁は手計算したことがありますが、 二度とやりたくありません 。

この記事を読む前に、 匿さんの某積分記事 を読んでおくとよいかもしれません。参考にはなりませんが。

円周率計算

概説

まず、次の等式は当然成り立ちます。

$$ \int_0^1\frac1{1+x^2}dx=\frac\pi4 $$

こんなの当たり前ですね。そこで、区間$(0,1)$においてかなり小さい値をもつ関数$f(x)$を用意して、

$$ \int_0^1\frac{f(x)}{1+x^2}dx,\int_0^1f(x)dx $$

の二つが高校範囲で計算できてかつ、表式に$\pi$と有理数のみが出てこればこちらのものです。具体的にはこんな不等式評価が成り立ちます。

$$ \frac12\int_0^1f(x)dx<\int_0^1\frac{f(x)}{1+x^2}dx<\int_0^1f(x)dx $$

証明はほぼ当たり前なので少し天下り的になりますが、結果的に次のような取り方をしました。

$$ f(x)=x^{4n}(1-x)^{4n} $$

はい。ベータ関数です。計算自体は(部分積分を繰り返し用いるなどすれば)可能なのですが、面倒なのでズルをします。

\begin{aligned} \int_0^1x^{4n}(1-x)^{4n}dx &=B(4n+1,4n+1)\\ &=\frac{(4n)!^2}{(8n+1)!}\\ &=\frac1{(8n+1)\binom{8n}{4n}} \end{aligned}

みていただくとわかる通り、分母に中心二項係数が出てくるので$n$が大きくなるとかなり急速に$0$に近づきます。細かいことが知りたければStirling近似とかを使えばいいですね。中の積分も面倒ではありますが、$\pi$と有理数しか出てこないことは技術的には簡単に示せるので、実質的問題は

$$ \int_0^1\frac{x^{4n}(1-x)^{4n}}{1+x^2}dx $$

を求める作業に帰着されました。

$n=1$の場合

参考程度に$n$が小さい場合について具体的に計算します。まず、ベータ積分のほうは

$$ B(5,5)=\frac{4!^2}{9!}=\frac1{630} $$

と計算できます。真ん中の積分は普通に割り算して計算してあげると

\begin{aligned} \int_0^1\frac{x^{4}(1-x)^{4}}{1+x^2}dx &=\int_0^1\left(x^6-4x^5+5x^4-4x^2+4-\frac4{1+x^2}\right)\\ &=\left(\frac17-\frac46+1-\frac43+4\right)-\pi\\ &=\frac{22}7-\pi \end{aligned}

であるので、次の不等式を得ます。

$$ \frac1{1260}<\frac{22}7-\pi<\frac1{630} $$

これを整理すると、

$$ \frac{22}7-\frac1{630}<\pi<\frac{22}7-\frac1{1260} $$

となり、この結果だけで余裕で小学校における定説(?)である$\pi=3.14$を算出できます。この計算はさほど大変ではありません。恐らく、手計算で実用的なレベルの評価としては$n=2$くらいが限界な気がします。

$n=2$の有理数評価もしておくと、地道な計算により
\begin{aligned}B(9,9)&=\frac1{218790},\\\int_0^1\frac{x^{8}(1-x)^{8}}{1+x^2}dx&=4\pi-\frac{188684}{15015}\end{aligned}
なので、

$$ \frac1{1750320}+\frac{188684}{60060}<\pi<\frac1{875160}+\frac{188684}{60060} $$

であるから、頑張って割り算すると

$$ 3.14159231...<\pi<3.1415928... $$

によって有効数字7桁を決めることができます。これくらいの計算なら受験本番でもできていいと思います。ちなみに$\frac{355}{133}$よりも少しだけいい近似なくらいなので、覚えたいならむしろこっちのほうがいいですね。

本題

とりあえず精度が欲しいので、$n=25$くらいにして積分を計算してみます。これは一度だけ手計算したことがあります。当時の僕はこれだと余裕で$50$桁に届くことがStirling評価からわかるので、$50$桁だけ計算しました。この結果、なんと(高校一年生の時に教室の後ろの黒板に書いて)写メ取ってたみたいです。見にくいので有理数に関するところだけ赤くします。また、計算過程なんてものはやりたくもないので結果だけ書いていきます。まずベータ積分の方は、

$$ B(101,101)=\color{red}\frac1{18200251445876759514246239592574316938775422524758080705105320} $$

となります。それに対して、不等式の内側の積分は

$$ \int_0^1\frac{x^{100}(1-x)^{100}}{1+x^2}dx={\color{red}\frac{25422588845948287917840198990310955266893220891255517493881289524342496083962053135792429098585663}{28749487633387736328132736580036532695803192488438380846224993685496448911426147450}}-281474976710656\pi $$

にと(8時間くらい頑張って)計算できるので、結果的に不等式は次のような形になります。

\begin{aligned} &\quad\color{red}\frac1{36400502891753519028492479185148633877550845049516161410210640}\\ &<{\color{red}\frac{25422588845948287917840198990310955266893220891255517493881289524342496083962053135792429098585663}{28749487633387736328132736580036532695803192488438380846224993685496448911426147450}}-281474976710656\pi\\ &<\color{red}\frac1{18200251445876759514246239592574316938775422524758080705105320} \end{aligned}

これを両辺$281474976710656$で割って移項などの整理を行うと、（ここは赤くしてません）

\begin{aligned} &\frac{25422588845948287917840198990310955266893220891255517493881289524342496083962053135792429098585663}{8092261362051105765248285301685603166696861085875287087233705766880133617989065376859671442227200}-\frac1{10245830703712488149502452713532562564917398605977053684060123252079292579840}\\ &<\pi<\\ &\frac{25422588845948287917840198990310955266893220891255517493881289524342496083962053135792429098585663}{8092261362051105765248285301685603166696861085875287087233705766880133617989065376859671442227200}-\frac1{5122915351856244074751226356766281282458699302988526842030061626039646289920} \end{aligned}

となり、とりあえず凄そうな分数が出てくるみたいです。$50$桁を超える十進小数だけ何故か過去の研究ノートに結果が書いてないのでWolframに入れようかと思ったのですが...なんと、文字列長の制限を食らいました。式が長すぎるようです。そこで、こんな計算やってられないので こちらの計算機 を使うことにしました。ここからは手計算じゃないですよ。先の分数を十進小数に直した不等式にしてあげると、

\begin{aligned} &3.141592653589793238462643383279502884197169399375105820974944592307816406286267527922991584361811747032194196784852017799295308862\\ &<\pi<\\ &3.141592653589793238462643383279502884197169399375105820974944592307816406286169927247271102405380017296915582316302269556264686866 \end{aligned}

となるらしいので、最右辺と最左辺がどこまで一致しているか色を付けてみると、

\begin{aligned} \mathrm{LHS}&={\color{red}3.141592653589793238462643383279502884197169399375105820974944592307816406286}267527922991584361811747032194196784852017799295308862\\ \mathrm{RHS}&={\color{red}3.141592653589793238462643383279502884197169399375105820974944592307816406286}169927247271102405380017296915582316302269556264686866 \end{aligned}

ええと、これは何桁なんだ。数えるのも面倒だったので こちらのサイト に数えてもらったところ、$\color{red}76$桁の精度であるようだ。精度高くない？

まとめ

円周率の近似値は、おおよそ$3.141592653589793238462643383279502884197169399375105820974944592307816406286$であることがわかりました！

$\ln2$も計算してみる

$2$の自然対数も求めることができます。具体的に言うと

$$ \frac12\int_0^1x^n(1-x)^ndx<\int_0^1\frac{x^n(1-x)^n}{1+x}dx<\int_0^1x^n(1-x)^ndx $$

で不等式を作ればよいです。対数の引数の値が大きくなると挟む積分が極限で消えてくれないというバグが発生するので$\ln 5$などを求めるのには使えません。(個人的には、$2^{10}=1024$と$2^{13}=8192$を繰り返しかけて$10$のべき乗の近似値を作って$\operatorname{lc}2$を求めてから、間接的に$\ln5$を求めるほうが現実的だと思います)ここでは$n=1,2,3,4,8$のときの$\ln2$の有理数近似を計算します。まず積分の結果については

\begin{aligned} \int_0^1\frac{x(1-x)}{1+x}dx &=\frac32-2\ln2\\ \int_0^1\frac{x^2(1-x)^2}{1+x}dx &=4\ln2-\frac{11}4\\ \int_0^1\frac{x^3(1-x)^3}{1+x}dx&=\frac{111}{20}-8\ln2\\ \int_0^1\frac{x^4(1-x)^4}{1+x}dx&=16\ln2-\frac{621}{56}\\ \int_0^1\frac{x^8(1-x)^8}{1+x}dx&=256\ln2-\frac{42629549}{240240} \end{aligned}

になります。ベータ関数のほうは

\begin{aligned} B(2,2)&=\frac16\\ B(3,3)&=\frac1{30}\\ B(4,4)&=\frac1{140}\\ B(5,5)&=\frac1{630}\\ B(9,9)&=\frac1{218790} \end{aligned}

であることがわかるので、$n=1,2,3,4,8$の順に

\begin{aligned} \frac23&<\ln2<\frac34-\frac1{24}\\ \frac{11}{16}+\frac1{240}&<\ln2<\frac{11}{16}+\frac1{120}\\ \frac{111}{160}-\frac1{1120}&<\ln2<\frac{111}{160}-\frac1{2240}\\ \frac{621}{896}+\frac1{20160}&<\ln2<\frac{621}{896}+\frac1{10080}\\ \frac{42629549}{61501440}+\frac1{112020480}&<\ln2<\frac{42629549}{61501440}+\frac1{56010240} \end{aligned}

で、小数表示では大まかには(ノーテーションも雑になっているのは許してください)

\begin{aligned} 0.67&<\ln2<0.70\\ 0.6916&<\ln2<0.6958\\ 0.6928&<\ln2<0.6933\\ 0.69312&<\ln2<0.69317\\ 0.6931471775&<\ln2<0.6931471864 \end{aligned}

なので、$\ln2$の近似値が$0.69314717$であることが分かりました。$n=4$のときによく知られている$0.6931$まで計算できるのは初めて計算した時は驚きましたね。

おわりに

一番最初にこの記事を書こうと思ったのは この記事 までさかのぼりますが、研究ノートの発掘が面倒だったのでずっとさぼっていました。どうせツイートしてしまったので、作ってみるか、という動機で(対数のほうの)計算をしていたのですが、$n=8$に関してはさすがにやったことがなく、かなり高い精度で計算できたというのが少し面白かったです。

近似計算が数学において本質的に重要か、と言われれば必ずしもそうでないと思いますが、何処までの精度をもって計算できるかというのはときたま意外な一面を見せるので面白いですね。あと誰か、この方法(あるいはより一般的なもの)を知っている人がいたら教えてください。

最後まで読んでいただきありがとうございました。