3
大学数学基礎解説
文献あり

円周率を10進数で任意の桁数の精度で近似できる式で遊ぼう(日曜数学会)

1078
0

はじめに

AIが描いた 円周率を計算する少女たちの図(使用ツール:イラストお絵描きばりぐっどくん) AIが描いた 円周率を計算する少女たちの図(使用ツール:イラストお絵描きばりぐっどくん)

 この記事は、2023211 日にオンラインで開催された 第26回日曜数学会 で私が発表した内容に大幅に加筆修正したものになります。

・ 発表時のスライド資料はこちら→ 【スライド資料】

概要

 先月(20231 月)Twitterで話題になった「円周率の10進数n桁目を直接計算できる式」は、実は伝言ゲームの結果、遊び方を誤って拡散されてしまった式です。この記事では、私の考える「この式の正しい遊び方」をご紹介します。

円周率のn桁目が得られる式

 まず、次の式をご覧ください。

【Twitterで拡散された式】
  π(2(1)n+1(2n)!22nB2n(122n)(132n)(152n)(172n))1/(2n)

 (拡散の過程での誤植部分は修正しています。)

 この式はTwitterで「10進数での n 桁目を直接(それまでの桁を計算せずに)計算できる式」という触れ込みで拡散されました。
 (ここからは、この式のことを単に「例の式」と書きます。)

 例の式の元ネタとなった論文を確認したところ、どうやらそれは事実ではないことに気が付きました。

【元ネタとなった論文: Simon Plouffe,A formula for the n-th decimal digit or binary of π and powers of π,2022

 実際には例の式は「10進数で 2n 桁の精度で円周率を近似する式」でした。『直接(それまでの桁を計算せずに)』計算することはできませんが、欲しい精度の円周率の近似値がラクに得られるので、任意の桁の数字を『効率よく』得られる式となります。

 もう少し詳しく説明しましょう。

 この式の中にある** B2n **は「関・ベルヌーイ数」と呼ばれる有理数です。

 【偶数番目の関・ベルヌーイ数の表】

B2B4B6B8B10B12B14B16
161301421305666912730763617510

 したがって、例の式の両辺を 2n 乗すると、 π2n を近似する有理数が得られます!

【例の式の両辺を 2 乗した式】
  π2n2(1)n+1(2n)!22nB2n(122n)(132n)(152n)(172n)

 ・ 右辺は有理数になっているので計算が簡単!

 このとき、近似の精度は 10 進数でおよそ 2n 桁以上の精度があります!
 その理由は今から説明します。

解説(前半)

リーマンゼータ関数

 例の式を作るのに最重要な道具は、ズバリ、「リーマンゼータ関数」です!

 リーマンゼータ関数というのは、次のような式で定義される関数です。数論をはじめ数学や物理学など様々な分野で用いられている関数です。「見たことある!」という人も多いのではないでしょうか。  

リーマンゼータ関数

ζ(s):=n=11ns=11s+12s+13s+14s+

 リーマンゼータ関数は 𝑠 が大きくなるとどんどん 1 に近づいていきます。式の形を見れば、その理由はなんとなくわかると思います。

1 へと近づく

ζ(2)=112+122+132+142+=1.64493

ζ(4)=114+124+134+144+=1.08232

ζ(6)=116+126+136+146+=1.01734

ζ(8)=118+128+138+148+=1.00407

ζ(10)=1110+1210+1310+1410+=1.00099

ζ(12)=1112+1212+1312+1412+=1.00024

(注) この表で s が偶数の場合しか書いていないのは、単にこの記事では偶数の場合しか使わないからです。

 また、s が偶数のときは、その厳密値が次の表のような「有理数 ×πs」の形で得られることも知られています。これらの値を「特殊値」と呼ぶこともあります。

s が正の偶数の場合の特殊値

ζ(2)=π26

ζ(4)=π490

ζ(6)=π6945

ζ(8)=π89450

ζ(10)=π1093555

ζ(12)=691π12638512875

円周率の偶数乗

さて、上の表(例 2)の式を少し変形するとこのような形にできます。

πs をリーマンゼータ関数で表す

π2=6ζ(2)

π4=90ζ(4)

π6=945ζ(6)

π8=9450ζ(8)

π10=93555ζ(10)

π12=638512875691ζ(12)

 s が大きいときは ζ(s)1 となることから、この表の ζ(s)1 に置き換えても、πs のざっくりとした近似にはなります。が、このままでは精度はそれほどよくありません。

リーマンゼータ関数のオイラー積表示

 そこで、ζ(s) を簡単な有理数で近似します。リーマンゼータ関数は、「オイラー積表示」と言って、素数を使った無限積の形で表すことができるのですが、これを使います。

リーマンゼータ関数のオイラー積表示

ζ(s)=p;prime11ps=112s113s115s117s1111s

 ただし、右辺の積記号はすべての素数の積を表す。

導出方法

クリック/タップで詳細
 積表示は、素因数分解の一意性を使って因数分解の形に変形することで導出できます。

ζ(s)=11s+12s+13s+14s+15s+16s+=11s+12s+13s+1(22)s+15s+1(23)s+=(1+12s+122s+123s+)×(1+13s+132s+133s+)×(1+15s+152s+153s+)=112s113s115s117s1111s(等比級数の和の公式を使った)=p;prime11ps

 くわしくは、tsujimotterさんのサイトをご覧ください。

  tsujimotterのノートブック:ゼータ関数のオイラー積

 次に、オイラー積表示を第 4 項までで打ち切ることで近似した式を考えましょう。

リーマンゼータ関数のオイラー積表示を第 4 項までで打ち切った近似式

ζ(s)112s113s115s117s

 これを使って具体的に πs の近似値を計算するとこうなります。

πs を近似する有理数

π26(122)(132)(152)(172)=91.57031

π490(124)(134)(154)(174)=972.39633

π6945(126)(136)(156)(176)=961.385838

π89450(128)(138)(158)(178)=9488.536095

π1093555(1210)(1310)(1510)(1710)=93648.04747101

π12638512875691(1212)(1312)(1512)(1712)=924269.181523120

 右辺の値の下に、実際の値と何桁まで一致しているか表示してあります。
 おおよそ、πs10 進数 s 桁くらいの精度で近似できていることがわかります。
 オイラー積の第 5 項以後、すなわち、11 以上の素数に関する項を省略しても大丈夫なのはなぜでしょうか。
 考えてみると、11 以上の素数に関する項は、計算結果に 10s 以下の影響しか与えません。実際、11 以上の素数に関する項を全て省略したとしても計算結果に与える影響は 10s より小さいと考えられます。(厳密には証明が必要ですが、ここでは数値計算の結果で判断しています。)

1/(2n) 乗すると誤差が小さくなる

 つぎに、上の表(例 5)の両辺を 1/(2n) 乗することで π の近似値を計算してみましょう。

π(6(122)(132)(152)(172))1/2=31.09359

π(90(124)(134)(154)(174))1/4=3.141448

π(945(126)(136)(156)(176))1/6=3.14159272

π(9450(128)(138)(158)(178))1/8=3.1415926591

π(93555(1210)(1310)(1510)(1710))1/10=3.1415926535117

π(638512875691(1212)(1312)(1512)(1712))1/12=3.141592653589136

 さっきより精度があがりました。ちょっと見にくいかもしれませんので、表にしてみましょう。

nπ2n の近似値の精度π の近似値の精度
111
224
357
469
51011
61213

 1/(2n) 乗することで精度があがり、特に、n2 以上のときは、π2n 桁以上の精度で求めることができていることがわかりますね。

なぜ精度が上がるのか

 n を自然数、a1 より大きい実数 とし、a に近い実数 x を考えます。 xnan の差が Δ とします。つまり

 xn=an+Δ

 このとき、n が十分大きければ、xa の差(の絶対値)は、 Δ (の絶対値)よりも必ず小さくなります。このことは、

 xa+Δnan1

と近似できることからわかります。

n=50 の場合

 ちなみに、n=50 だと何桁くらいの精度になるでしょうか?
 実際に調べてみましょう。


ζ(100)=189196075638244250590454866138987443745405683066133872266771392408622790830394495422π1009815205420757514710108178059369553458327392260750404049930407987933582359080767225644716670683512153512547802166033089160919189453125

ですから、


π(9815205420757514710108178059369553458327392260750404049930407987933582359080767225644716670683512153512547802166033089160919189453125189196075638244250590454866138987443745405683066133872266771392408622790830394495422(12100)(13100)(15100)(17100))1/100=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982141057

となり、2n(=100) 桁を上回る、105 桁の精度で π の近似値を計算することができました。

 例の式が 10 進数 2n 桁の精度で円周率が求められる理由は、これで説明できました。

解説(後半)

 ところで、例の式の中にある「関・ベルヌーイ数」についてまだ説明していませんでした。
 というわけで、ここからは関・ベルヌーイ数について説明していきます。

関・ベルヌーイ数

 関・ベルヌーイ数というのは、べき乗和の公式と深い関係のある有理数の数列です。(余談ですが、「数列」なのに「~数」と呼ぶのはなぜでしょうね?)

 具体的には、次のような関係になります。

10+20+30++n0=11B0n111+21+31++n1=12((20)B0n2+(21)B1n1)12+22+32++n2=13((30)B0n3+(31)B1n2+(32)B2n1)j=1njk=1k+1j=0k(k+1j)Bjnkj+1

(注) 関・ベルヌーイ数については 2 とおりの定義がありますが、この記事では深入りしません。

リーマンゼータ関数と関・ベルヌーイ数の関係

 さきほど、リーマンゼータ関数の引数が正の偶数のときは、「有理数 ×πs」の形で書くことができる、と書きました。実は、この有理数の部分を、関・ベルヌーイ数を使って表すことができることがしられています。
 具体的には、次の様に表すことができます。

リーマンゼータ関数の特殊値を関・ベルヌーイ数を使って表す

任意の正の偶数 2n に対して、

ζ(2n)=(1)n+1B2n(2π)2n2(2n)!

 この式の証明は簡単ではありませんが、例えば、次のサイトなどを参考にしてください。

青空学園数学科書庫:ζ(2l)を関・ベルヌーイ数で表す

 これを使って先ほどの式を書き換えることで例の式が得られます!

例の式(再掲)

  π(2(1)n+1(2n)!22nB2n(122n)(132n)(152n)(172n))1/(2n)

応用 平方根で計算できるようにする

 例の式の仕組みがわかったことで、ここからは応用を考えます。

 こんなことありませんか?「急に円周率を計算してみたくなったなー。精度はそんなになくても、簡単に計算できる方法ないかな?

 実は、例の式はもう少し簡単にできます!
 1/(2n) 乗を 1/2 乗(つまり平方根)にして、階乗や2のべき乗を使わない形にできます!

 つまり、少しくらいの桁なら、がんばれば手計算できるくらいの簡単な形にすることができます!

(「手計算できる」とはいえ、ベルヌーイ数はあらかじめ知っておく必要はあります💦)

2 つの式から新しい式を作る

例の式の nn+1 にしたものと、元の式を並べます。

π2n+2=2(1)n(2n+2)!22n+2B2n+2(122n2)(132n2)(152n2)(172n2)

π2n=2(1)n+1(2n)!22nB2n(122n)(132n)(152n)(172n)

 2 つの式の辺々を割り算して整理すると

π2B2n(2n+2)(2n+1)(22n+222)(32n+232)(52n+252)(72n+272)4B2n+2(22n+21)(32n+21)(52n+21)(72n+21)

 両辺の平方根をとって

πB2n(2n+2)(2n+1)(22n+222)(32n+232)(52n+252)(72n+272)4B2n+2(22n+21)(32n+21)(52n+21)(72n+21)

改良版 π の近似値

πB2n(2n+2)(2n+1)(22n+222)(32n+232)(52n+252)(72n+272)4B2n+2(22n+21)(32n+21)(52n+21)(72n+21)

 1/(2n) 乗の部分が平方根になりました。さらに、2 のべき乗や階乗などが消えて、計算しやすくなりました!
 平方根は筆算で開く方法が知られていますから、少しくらいの桁ならがんばれば手計算できますね!
 
・  高校数学の美しい物語:開平法のやり方と原理

・  手計算の例

n=1,2,3,50 のとき

B243(2422)(3432)(5452)(7472)4B4(241)(341)(541)(741)=3.129

B465(2622)(3632)(5652)(7672)4B6(261)(361)(561)(761)=3.14147

B687(2822)(3832)(5852)(7872)4B8(281)(381)(581)(781)=3.1415963


B100102101(210222)(310232)(510252)(710272)4B102(21021)(31021)(51021)(71021)=3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798211045

 ちょうど 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 )

改良版の式と元論文との関係

 上記の改良版の式は、元ネタとなった論文にはない式です。
 ……が、全くのオリジナルというわけではありません(告白)。

 実は、論文にはこんな式がでていました。

【論文にあった式】

π212B2n(n+1)(2n+1)B2n+2

(注) 記事の都合上符号を修正しています。

 改良版の式とよく似ていますが、リーマンゼータ関数の部分を 省略している点が異なります。
 このままでは精度があまりよくないので、リーマンゼータ関数の近似式を付け加えて改良版の式をつくったのでした。  

おわりに

 というわけで、例の式の私の見つけた遊び方を紹介してみました。
 拡散のされ方でちょっとケチがついてしまいましたが、「10進数で n 桁くらいの精度で円周率が計算したいな~」というニッチな願いが叶えられるステキな式だと思います。
 ほかにもいろいろな遊び方ができると思いますので皆さんも遊んでみてください!

【例の式】

  π(2(1)n+1(2n)!22nB2n(122n)(132n)(152n)(172n))1/(2n)

【改良版】

πB2n(2n+2)(2n+1)(22n+222)(32n+232)(52n+252)(72n+272)4B2n+2(22n+21)(32n+21)(52n+21)(72n+21)

【関・ベルヌーイ数】

B2B4B6B8B10B12B14B16
161301421305666912730763617510

参考文献

投稿日:2023215
OptHub AI Competition

この記事を高評価した人

高評価したユーザはいません

この記事に送られたバッジ

バッジはありません。
バッチを贈って投稿者を応援しよう

バッチを贈ると投稿者に現金やAmazonのギフトカードが還元されます。

投稿者

apu_yokai
apu_yokai
485
65578

コメント

他の人のコメント

コメントはありません。
読み込み中...
読み込み中
  1. はじめに
  2. 概要
  3. 円周率のn桁目が得られる式
  4. 解説(前半)
  5. リーマンゼータ関数
  6. 円周率の偶数乗
  7. リーマンゼータ関数のオイラー積表示
  8. 1/(2n) 乗すると誤差が小さくなる
  9. なぜ精度が上がるのか
  10. n=50 の場合
  11. 解説(後半)
  12. 関・ベルヌーイ数
  13. リーマンゼータ関数と関・ベルヌーイ数の関係
  14. 応用 平方根で計算できるようにする
  15. 2 つの式から新しい式を作る
  16. 改良版の式と元論文との関係
  17. おわりに
  18. 参考文献