AIが描いた 円周率を計算する少女たちの図(使用ツール:イラストお絵描きばりぐっどくん)
この記事は、$2023$年 $2$ 月 $11$ 日にオンラインで開催された 第26回日曜数学会 で私が発表した内容に大幅に加筆修正したものになります。
・ 発表時のスライド資料はこちら→ 【スライド資料】
先月($2023$ 年 $1$ 月)Twitterで話題になった「円周率の10進数n桁目を直接計算できる式」は、実は伝言ゲームの結果、遊び方を誤って拡散されてしまった式です。この記事では、私の考える「この式の正しい遊び方」をご紹介します。
まず、次の式をご覧ください。
【Twitterで拡散された式】
${\displaystyle
\pi\simeq\left(
\frac{2(-1)^{n+1}(2n)!}
{2^{2n}B_{2n}(1-2^{-2n})(1-3^{-2n})(1-5^{-2n})(1-7^{-2n})}
\right)^{1/(2n)}
}$
(拡散の過程での誤植部分は修正しています。)
この式はTwitterで「10進数での n 桁目を直接(それまでの桁を計算せずに)計算できる式」という触れ込みで拡散されました。
(ここからは、この式のことを単に「例の式」と書きます。)
例の式の元ネタとなった論文を確認したところ、どうやらそれは事実ではないことに気が付きました。
【元ネタとなった論文: Simon Plouffe,A formula for the n-th decimal digit or binary of π and powers of π,2022 】
実際には例の式は「10進数で 2n 桁の精度で円周率を近似する式」でした。『直接(それまでの桁を計算せずに)』計算することはできませんが、欲しい精度の円周率の近似値がラクに得られるので、任意の桁の数字を『効率よく』得られる式となります。
もう少し詳しく説明しましょう。
この式の中にある** $B_{2n}$ **は「関・ベルヌーイ数」と呼ばれる有理数です。
【偶数番目の関・ベルヌーイ数の表】
$B_2$ | $B_4$ | $B_6$ | $B_8$ | $B_{10}$ | $B_{12}$ | $B_{14}$ | $B_{16}$ | $\cdots$ |
---|---|---|---|---|---|---|---|---|
$\frac{1}{6}$ | $\frac{-1}{30}$ | $\frac{1}{42}$ | $\frac{-1}{30}$ | $\frac{5}{66}$ | $\frac{-691}{2730}$ | $\frac{7}{6}$ | $\frac{-3617}{510}$ | $\cdots$ |
したがって、例の式の両辺を $2n$ 乗すると、 $\pi^{2n}$ を近似する有理数が得られます!
【例の式の両辺を $2$ 乗した式】
${\displaystyle
\pi^{2n}\simeq
\frac{2(-1)^{n+1}(2n)!}
{2^{2n}B_{2n}(1-2^{-2n})(1-3^{-2n})(1-5^{-2n})(1-7^{-2n})}
}$
・ 右辺は有理数になっているので計算が簡単!
このとき、近似の精度は $10$ 進数でおよそ $2n$ 桁以上の精度があります!
その理由は今から説明します。
例の式を作るのに最重要な道具は、ズバリ、「リーマンゼータ関数」です!
リーマンゼータ関数というのは、次のような式で定義される関数です。数論をはじめ数学や物理学など様々な分野で用いられている関数です。「見たことある!」という人も多いのではないでしょうか。
${\displaystyle \zeta (s):=\sum _{n=1}^{\infty }{\dfrac {1}{n^{s}}}={\dfrac {1}{1^{s}}}+{\dfrac {1}{2^{s}}}+{\dfrac {1}{3^{s}}}+{\dfrac {1}{4^{s}}}+\cdots }$
リーマンゼータ関数は 𝑠 が大きくなるとどんどん $1$ に近づいていきます。式の形を見れば、その理由はなんとなくわかると思います。
$\zeta(2)=\dfrac{1}{1^2}+\dfrac{1}{2^2}+\dfrac{1}{3^2}+\dfrac{1}{4^2}+\cdots=1.64493\cdots$
$\zeta(4)=\dfrac{1}{1^4}+\dfrac{1}{2^4}+\dfrac{1}{3^4}+\dfrac{1}{4^4}+\cdots=1.08232\cdots$
$\zeta(6)=\dfrac{1}{1^6}+\dfrac{1}{2^6}+\dfrac{1}{3^6}+\dfrac{1}{4^6}+\cdots=1.01734\cdots$
$\zeta(8)=\dfrac{1}{1^8}+\dfrac{1}{2^8}+\dfrac{1}{3^8}+\dfrac{1}{4^8}+\cdots=1.00407\cdots$
$\zeta(10)=\dfrac{1}{1^{10}}+\dfrac{1}{2^{10}}+\dfrac{1}{3^{10}}+\dfrac{1}{4^{10}}+\cdots=1.00099\cdots$
$\zeta(12)=\dfrac{1}{1^{12}}+\dfrac{1}{2^{12}}+\dfrac{1}{3^{12}}+\dfrac{1}{4^{12}}+\cdots=1.00024\cdots$
(注) この表で $s$ が偶数の場合しか書いていないのは、単にこの記事では偶数の場合しか使わないからです。
また、$s$ が偶数のときは、その厳密値が次の表のような「有理数 $\times \pi^{s}$」の形で得られることも知られています。これらの値を「特殊値」と呼ぶこともあります。
$\zeta(2)=\dfrac{\pi^{2}}{6}$
$\zeta(4)=\dfrac{\pi^{4}}{90}$
$\zeta(6)=\dfrac{\pi^{6}}{945}$
$\zeta(8)=\dfrac{\pi^{8}}{9450}$
$\zeta(10)=\dfrac{\pi^{10}}{93555}$
$\zeta(12)=\dfrac{691\pi^{12}}{638512875}$
さて、上の表(例 $2$)の式を少し変形するとこのような形にできます。
$\pi^{2}=6\,\zeta(2)$
$\pi^{4}=90\,\zeta(4)$
$\pi^{6}=945\,\zeta(6)$
$\pi^{8}=9450\,\zeta(8)$
$\pi^{10}=93555\,\zeta(10)$
$\pi^{12}=\dfrac{638512875}{691}\,\zeta(12)$
$s$ が大きいときは $\zeta(s)\simeq 1$ となることから、この表の $\zeta(s)$ を $1$ に置き換えても、$\pi^s$ のざっくりとした近似にはなります。が、このままでは精度はそれほどよくありません。
そこで、$\zeta(s)$ を簡単な有理数で近似します。リーマンゼータ関数は、「オイラー積表示」と言って、素数を使った無限積の形で表すことができるのですが、これを使います。
${\displaystyle \begin{align} \zeta(s) &=\prod_{p;prime}\frac{1}{1-p^{-s}}\\ &=\frac{1}{1-2^{-s}}\cdot\frac{1}{1-3^{-s}}\cdot\frac{1}{1-5^{-s}}\cdot\frac{1}{1-7^{-s}}\cdot\frac{1}{1-11^{-s}}\cdots\\ \end{align} }$
ただし、右辺の積記号はすべての素数の積を表す。
次に、オイラー積表示を第 $4$ 項までで打ち切ることで近似した式を考えましょう。
${\displaystyle \begin{align} \zeta(s) &\simeq\frac{1}{1-2^{-s}}\cdot\frac{1}{1-3^{-s}}\cdot\frac{1}{1-5^{-s}}\cdot\frac{1}{1-7^{-s}} \end{align} }$
これを使って具体的に $\pi^s$ の近似値を計算するとこうなります。
$\pi^{2}\simeq\dfrac{6}{(1-2^{-2})(1-3^{-2})(1-5^{-2})(1-7^{-2})}=\underbrace{9}_{1桁}.57031\cdots$
$\pi^{4}\simeq\dfrac{90}{(1-2^{-4})(1-3^{-4})(1-5^{-4})(1-7^{-4})}=\underbrace{97}_{2桁}.39633\cdots$
$\pi^{6}\simeq\dfrac{945}{(1-2^{-6})(1-3^{-6})(1-5^{-6})(1-7^{-6})}=\underbrace{961.38}_{5桁}838\cdots$
$\pi^{8}\simeq\dfrac{9450}{(1-2^{-8})(1-3^{-8})(1-5^{-8})(1-7^{-8})}=\underbrace{9488.53}_{6桁}095\cdots$
$\pi^{10}\simeq\dfrac{93555}{(1-2^{-10})(1-3^{-10})(1-5^{-10})(1-7^{-10})}=\underbrace{93648.04747}_{10桁}1\cdots$
$\pi^{12}\simeq\dfrac{638512875}{691(1-2^{-12})(1-3^{-12})(1-5^{-12})(1-7^{-12})}=\underbrace{924269.181523}_{12桁}0\cdots$
右辺の値の下に、実際の値と何桁まで一致しているか表示してあります。
おおよそ、$\pi^s$ を $10$ 進数 $s$ 桁くらいの精度で近似できていることがわかります。
オイラー積の第 $5$ 項以後、すなわち、$11$ 以上の素数に関する項を省略しても大丈夫なのはなぜでしょうか。
考えてみると、$11$ 以上の素数に関する項は、計算結果に $10^{-s}$ 以下の影響しか与えません。実際、$11$ 以上の素数に関する項を全て省略したとしても計算結果に与える影響は $10^{-s}$ より小さいと考えられます。(厳密には証明が必要ですが、ここでは数値計算の結果で判断しています。)
つぎに、上の表(例 $5$)の両辺を $1/(2n)$ 乗することで $\pi$ の近似値を計算してみましょう。
$\pi\simeq\left(\dfrac{6}{(1-2^{-2})(1-3^{-2})(1-5^{-2})(1-7^{-2})}\right)^{1/2}=\underbrace{3}_{1桁}.09359\cdots$
$\pi\simeq\left(\dfrac{90}{(1-2^{-4})(1-3^{-4})(1-5^{-4})(1-7^{-4})}\right)^{1/4}=\underbrace{3.141}_{4桁}48\cdots$
$\pi\simeq\left(\dfrac{945}{(1-2^{-6})(1-3^{-6})(1-5^{-6})(1-7^{-6})}\right)^{1/6}=\underbrace{3.141592}_{7桁}2\cdots$
$\pi\simeq\left(\dfrac{9450}{(1-2^{-8})(1-3^{-8})(1-5^{-8})(1-7^{-8})}\right)^{1/8}=\underbrace{3.14159265}_{9桁}1\cdots$
$\pi\simeq\left(\dfrac{93555}{(1-2^{-10})(1-3^{-10})(1-5^{-10})(1-7^{-10})}\right)^{1/10}=\underbrace{3.1415926535}_{11桁}7\cdots$
$\pi\simeq\left(\dfrac{638512875}{691(1-2^{-12})(1-3^{-12})(1-5^{-12})(1-7^{-12})}\right)^{1/12}=\underbrace{3.141592653589}_{13桁}6\cdots$
さっきより精度があがりました。ちょっと見にくいかもしれませんので、表にしてみましょう。
$n$ | $\pi^{2n}$ の近似値の精度 | $\pi$ の近似値の精度 |
---|---|---|
$1$ | $1$桁 | $1$桁 |
$2$ | $2$桁 | $4$桁 |
$3$ | $5$桁 | $7$桁 |
$4$ | $6$桁 | $9$桁 |
$5$ | $10$桁 | $11$桁 |
$6$ | $12$桁 | $13$桁 |
$1/(2n)$ 乗することで精度があがり、特に、$n$ が $2$ 以上のときは、$\pi$ を $2n$ 桁以上の精度で求めることができていることがわかりますね。
$n$ を自然数、$a$を $1$ より大きい実数 とし、$a$ に近い実数 $x$ を考えます。 $x^n$ と $a^n$ の差が $\Delta$ とします。つまり
$x^{n}=a^{n}+\Delta$
このとき、$n$ が十分大きければ、$x$ と $a$ の差(の絶対値)は、 $\Delta$ (の絶対値)よりも必ず小さくなります。このことは、
$x\simeq a + \frac{\Delta}{na^{n-1}}$
と近似できることからわかります。
ちなみに、$n=50$ だと何桁くらいの精度になるでしょうか?
実際に調べてみましょう。
ですから、
となり、$2n(=100)$ 桁を上回る、$105$ 桁の精度で $\pi$ の近似値を計算することができました。
例の式が $10$ 進数 $2n$ 桁の精度で円周率が求められる理由は、これで説明できました。
ところで、例の式の中にある「関・ベルヌーイ数」についてまだ説明していませんでした。
というわけで、ここからは関・ベルヌーイ数について説明していきます。
関・ベルヌーイ数というのは、べき乗和の公式と深い関係のある有理数の数列です。(余談ですが、「数列」なのに「~数」と呼ぶのはなぜでしょうね?)
具体的には、次のような関係になります。
${\displaystyle \begin{align} 1^0+2^0+3^0+\cdots+n^0 &=\frac{1}{1} B_0 n^1\\ 1^1+2^1+3^1+\cdots+n^1 &=\frac{1}{2}\left( {2 \choose 0}B_0 n^2 +{2 \choose 1}B_1 n^1 \right)\\ 1^2+2^2+3^2+\cdots+n^2 &=\frac{1}{3}\left( {3 \choose 0}B_0 n^3 +{3 \choose 1}B_1 n^2 +{3 \choose 2}B_2 n^1 \right)\\ \qquad\vdots&\qquad\vdots\\ \sum_{j=1}^{n} j^k &=\frac{1}{k+1}\sum_{j=0}^{k} {k+1 \choose j}B_j n^{k-j+1} \\ \end{align} }$
(注) 関・ベルヌーイ数については $2$ とおりの定義がありますが、この記事では深入りしません。
さきほど、リーマンゼータ関数の引数が正の偶数のときは、「有理数 $\times \pi^s$」の形で書くことができる、と書きました。実は、この有理数の部分を、関・ベルヌーイ数を使って表すことができることがしられています。
具体的には、次の様に表すことができます。
任意の正の偶数 $2n$ に対して、
${\displaystyle \zeta (2n)=(-1)^{n+1}\,{\frac {\,B_{2n}(2\pi )^{2n}}{2(2n)!}}}$
この式の証明は簡単ではありませんが、例えば、次のサイトなどを参考にしてください。
これを使って先ほどの式を書き換えることで例の式が得られます!
${\displaystyle \pi\simeq\left( \frac{2(-1)^{n+1}(2n)!} {2^{2n}B_{2n}(1-2^{-2n})(1-3^{-2n})(1-5^{-2n})(1-7^{-2n})} \right)^{1/(2n)} }$
例の式の仕組みがわかったことで、ここからは応用を考えます。
こんなことありませんか?「急に円周率を計算してみたくなったなー。精度はそんなになくても、簡単に計算できる方法ないかな?」
実は、例の式はもう少し簡単にできます!
1/(2n) 乗を 1/2 乗(つまり平方根)にして、階乗や2のべき乗を使わない形にできます!
つまり、少しくらいの桁なら、がんばれば手計算できるくらいの簡単な形にすることができます!
(「手計算できる」とはいえ、ベルヌーイ数はあらかじめ知っておく必要はあります💦)
例の式の $n$ を $n+1$ にしたものと、元の式を並べます。
${\displaystyle \pi^{2n+2}= \frac{2(-1)^{n}(2n+2)!} {2^{2n+2}B_{2n+2}(1-2^{-2n-2})(1-3^{-2n-2})(1-5^{-2n-2})(1-7^{-2n-2})} }$
${\displaystyle \pi^{2n}= \frac{2(-1)^{n+1}(2n)!} {2^{2n}B_{2n}(1-2^{-2n})(1-3^{-2n})(1-5^{-2n})(1-7^{-2n})} }$
$2$ つの式の辺々を割り算して整理すると
${\displaystyle \begin{align} \pi^{2}&\simeq \frac{-B_{2n}(2n+2)(2n+1)(2^{2n+2}-2^{2})(3^{2n+2}-3^{2})(5^{2n+2}-5^{2})(7^{2n+2}-7^{2})} {4B_{2n+2}(2^{2n+2}-1)(3^{2n+2}-1)(5^{2n+2}-1)(7^{2n+2}-1)} \end{align} }$
両辺の平方根をとって
${\displaystyle \begin{align} \therefore \pi &\simeq\sqrt{ \frac{-B_{2n}(2n+2)(2n+1)(2^{2n+2}-2^{2})(3^{2n+2}-3^{2})(5^{2n+2}-5^{2})(7^{2n+2}-7^{2})} {4B_{2n+2}(2^{2n+2}-1)(3^{2n+2}-1)(5^{2n+2}-1)(7^{2n+2}-1)} } \end{align} }$
${\displaystyle \begin{align} \pi &\simeq\sqrt{ \frac{-B_{2n}(2n+2)(2n+1)(2^{2n+2}-2^{2})(3^{2n+2}-3^{2})(5^{2n+2}-5^{2})(7^{2n+2}-7^{2})} {4B_{2n+2}(2^{2n+2}-1)(3^{2n+2}-1)(5^{2n+2}-1)(7^{2n+2}-1)} } \end{align} }$
$1/(2n)$ 乗の部分が平方根になりました。さらに、$2$ のべき乗や階乗などが消えて、計算しやすくなりました!
平方根は筆算で開く方法が知られていますから、少しくらいの桁ならがんばれば手計算できますね!
・
高校数学の美しい物語:開平法のやり方と原理
・ 手計算の例
${\displaystyle \begin{align} &\sqrt{ \frac{-B_{2}\cdot 4 \cdot 3\cdot(2^{4}-2^{2})(3^{4}-3^{2})(5^{4}-5^{2})(7^{4}-7^{2})} {4B_{4}(2^{4}-1)(3^{4}-1)(5^{4}-1)(7^{4}-1)} }\\ &=\underbrace{3.1}_{2桁} 9\cdots \end{align} }$
${\displaystyle \begin{align} &\sqrt{ \frac{-B_{4}\cdot 6 \cdot 5\cdot(2^{6}-2^{2})(3^{6}-3^{2})(5^{6}-5^{2})(7^{6}-7^{2})} {4B_{6}(2^{6}-1)(3^{6}-1)(5^{6}-1)(7^{6}-1)} }\\ &=\underbrace{3.141}_{4桁} 7\cdots \end{align} }$
${\displaystyle \begin{align} &\sqrt{ \frac{-B_{6}\cdot 8 \cdot 7\cdot(2^{8}-2^{2})(3^{8}-3^{2})(5^{8}-5^{2})(7^{8}-7^{2})} {4B_{8}(2^{8}-1)(3^{8}-1)(5^{8}-1)(7^{8}-1)} }\\ &=\underbrace{3.14159}_{6桁} 3\cdots \end{align} }$
ちょうど $10$ 進数で $2n$ 桁くらいの精度になっていますね!
「手計算なんかやってられるか!」という人のために、Wolfram Alpha で計算するためのリンクもおいておきますね。
[ReplaceAll[((-BernoulliB[2n] (2n+2)(2n+1)(2^{2n+2}-4)(3^{2n+2}-9)(5^{2n+2}-25)(7^{2n+2}-49))/(4BernoulliB[2n+2] (2^{2n+2}-1)(3^{2n+2}-1)(5^{2n+2}-1)(7^{2n+2}-1)))^(1/2), {n -> 50}]]( https://ja.wolframalpha.com/input?i=ReplaceAll%5B%28%28-BernoulliB%5B2n%5D+%282n%2B2%29%282n%2B1%29%282%5E%7B2n%2B2%7D-4%29%283%5E%7B2n%2B2%7D-9%29%285%5E%7B2n%2B2%7D-25%29%287%5E%7B2n%2B2%7D-49%29%29%2F%284BernoulliB%5B2n%2B2%5D+%282%5E%7B2n%2B2%7D-1%29%283%5E%7B2n%2B2%7D-1%29%285%5E%7B2n%2B2%7D-1%29%287%5E%7B2n%2B2%7D-1%29%29%29%5E%281%2F2%29%2C+%7Bn+-%3E+50%7D%5D )
上記の改良版の式は、元ネタとなった論文にはない式です。
……が、全くのオリジナルというわけではありません(告白)。
実は、論文にはこんな式がでていました。
【論文にあった式】
${\displaystyle \begin{align} \pi^2 &\simeq-\frac{1}{2} \frac{B_{2n}(n+1)(2n+1)} {B_{2n+2}} \end{align} }$
(注) 記事の都合上符号を修正しています。
改良版の式とよく似ていますが、リーマンゼータ関数の部分を 省略している点が異なります。
このままでは精度があまりよくないので、リーマンゼータ関数の近似式を付け加えて改良版の式をつくったのでした。
というわけで、例の式の私の見つけた遊び方を紹介してみました。
拡散のされ方でちょっとケチがついてしまいましたが、「10進数で n 桁くらいの精度で円周率が計算したいな~」というニッチな願いが叶えられるステキな式だと思います。
ほかにもいろいろな遊び方ができると思いますので皆さんも遊んでみてください!
【例の式】
${\displaystyle \pi\simeq\left( \frac{2(-1)^{n+1}(2n)!} {2^{2n}B_{2n}(1-2^{-2n})(1-3^{-2n})(1-5^{-2n})(1-7^{-2n})} \right)^{1/(2n)} }$
【改良版】
${\displaystyle \begin{align} \pi &\simeq\sqrt{ \frac{-B_{2n}(2n+2)(2n+1)(2^{2n+2}-2^{2})(3^{2n+2}-3^{2})(5^{2n+2}-5^{2})(7^{2n+2}-7^{2})} {4B_{2n+2}(2^{2n+2}-1)(3^{2n+2}-1)(5^{2n+2}-1)(7^{2n+2}-1)} } \end{align} }$
【関・ベルヌーイ数】
$B_2$ | $B_4$ | $B_6$ | $B_8$ | $B_{10}$ | $B_{12}$ | $B_{14}$ | $B_{16}$ | $\cdots$ |
---|---|---|---|---|---|---|---|---|
$\frac{1}{6}$ | $\frac{-1}{30}$ | $\frac{1}{42}$ | $\frac{-1}{30}$ | $\frac{5}{66}$ | $\frac{-691}{2730}$ | $\frac{7}{6}$ | $\frac{-3617}{510}$ | $\cdots$ |