Go to Qiita Advent Calendar 2024 Top
LoginSignup
5
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@kyamaz

ChatGPTに円周率を計算するプログラムをお願いしてみた話

Last updated at Posted at 2024-04-07

こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz :frog: です。

本稿で活用したChatGPTは、ChatGPT-3.5 です。

はじめに

2022年のことですが、Googleの研究者が円周率を100兆桁目まで計算1しました。100兆桁目の数字は「0」だそうです。この情報を ChatGPT に答えて欲しくて、試しに聞いていました。残念ながら「現実的な時間内で計算することは不可能です。検索してね！」というような回答が返ってきました。
そこで ChatGPT は、

  • どのくらいの桁数の円周率を知っているのでしょうか。
  • 円周率を計算するアルゴリズムをどのくらい知っているのでしょうか。

という素朴な疑問が湧いてきました。この２点について、本稿で調べてみましょう。
次のリンクにあるように、同じような課題を調べられている記事や動画が幾つか見つかると思います。興味がある方は他の情報もご確認ください。

ChatGPT が答えた円周率は1000桁

１つめの疑問「ChatGPTはどのくらいの桁数の円周率を知っているか」については、シンプルに試せました。数えてみると1,000桁まで表示されました。

:frog:〔 円周率を可能な限り長く表示してください 〕

〔ChatGPTの回答〕
円周率は、無限の桁数を持つ無理数であり、通常は有限の桁数で表現されますが、いくつかの桁を表示することは可能です。円周率を可能な限り長く表示すると、以下のようになります。

表示された円周率

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198

ChatGPT が答えた円周率を求めるプログラム

２つめの疑問「ChatGPTは、円周率を計算するアルゴリズムをどのくらい知っているか」について、次の順番で質問して試しました。

:frog:〔 円周率を計算で求めるプログラムが知りたい 〕
:frog:〔 円周率を計算するための無限級数で収束速度が最も速いものを教えてください 〕
:frog:〔 ほかに最も速く収束する級数を教えてください 〕・・・★

※この最後の質問(★)を繰り返すと、異なるアルゴリズムを１つずつ教えてくれます。
※最初のプロンプトで円周率を計算するプログラムを聞いていますので、次のプロンプトに対してもプログラム付きで回答してくれました。但し、説明とプログラムが対応（アルゴリズムとプログラムが対応）できていないこともありました。

幾つかのアルゴリズムが提示されたあとは、これまでに回答された中にある回答をするようになります。つまり、この状態になるまで繰り返せば、ChatGPTが知っている”円周率を計算するためのアルゴリズム”のバリエーションが想定できると考えてもよいでしょう。

参考資料234を調べてみると、ChatGPTの回答に間違いが含まれておりましたが、それを正しながら整理してこのあと列挙していきます。

ライプニッツの公式

$\displaystyle \begin{equation} \frac{\pi}{4}=\frac{1}{1}-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\frac{1}{11}+\cdots \end{equation}$

これは、グレゴリー・ライプニッツ級数（逆正接関数 arctan x のマクローリン展開）である次式(A)に $x=1$ を代入したときの級数です。

$\displaystyle \begin{align} \arctan{x} &= \sum^{\infty}_{n=0} \frac{(-1)^n}{2n+1} x^{2n+1} = x - \frac{1}{3}x^{3} + \frac{1}{5}x^{5} - \frac{1}{7}x^{7} + \frac{1}{9}x^{9} - \cdots \tag{A} \\
\mathrm{arccot}{x} &= \frac{\pi}{2}-\arctan{x} \tag{B} \end{align}$

提示されたプログラム「ライプニッツの公式」 
import math                    # `math`ライブラリは必要ない

def gauss_arccot(precision):   # ガウスの公式ではなく、ライプニッツの公式の実装
    pi_over_2 = 0.0
    sign = 1
    for n in range(0, precision):
        pi_over_2 += sign / (2 * n + 1)
        sign = -sign
    return pi_over_2

def pi_from_gauss_arccot(precision):
    return 4 * gauss_arccot(precision)   # 出力されたプログラムは、"2*"でしたが"4*"が正しい

# 任意の桁数で円周率を計算
num_digits = int(input("計算したい円周率の桁数を入力してください: "))
precision = int(10 ** (num_digits / 2))
pi = pi_from_gauss_arccot(precision)
print("円周率:", pi)

マチンの公式

$\displaystyle \begin{equation} \frac{\pi}{4}= 4 \arctan{\frac{1}{5}} - \arctan{\frac{1}{239}} \end{equation}$

マチンの公式のWikiPediaによると、$\arctan{x}$の級数展開の式(A)を用いた方式は1970年代まで利用されていた方式だったようです。なお、2002年に日立のスーパーコンピューターを用いた式5も$\arctan{x}$の級数展開を使っていますが、ChatGPTはこの式を知らないようです。

提示されたプログラム「マチンの公式」

ChatGPTから「マチンの公式」の名称と説明はありましたが、プログラムはアルゴリズム通りの処理ではなく、対応するプログラムは提示されませんでした。

ベイリー＝ボールウェイン＝プラウフ（BBP）の公式

BBP（Bailey–Borwein–Plouffe formula）公式は次式で表されます。
$\displaystyle \pi = \sum^{\infty}_{k=0}(\frac{1}{16^{k}}(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}))$
このアルゴリズムは、先行する桁を計算せずに$\pi$の十六進法で $n$ 桁目（二進法で $4n$ 桁目）の数を直接求めることが特徴です。

提示されたプログラム「BBPアルゴリズム」 
import math

def bbp_algorithm(digits):
    pi = 0
    for k in range(digits):
        pi += 1 / 16 ** k * (4 / (8 * k + 1) - 2 / (8 * k + 4) - 1 / (8 * k + 5) - 1 / (8 * k + 6))
    return pi

# 求めたい円周率の桁数を入力
digits = int(input("計算したい円周率の桁数を入力してください: "))
pi = bbp_algorithm(digits)
print("円周率:", pi)

:frog:〔 mathライブラリではなく本来はdecimalライブラリを使って以下のように実装するのがよいです 〕

from decimal import *

def bbp_algorithm(digits):
    getcontext().prec = digits + 1
    pi = Decimal(0)
    for k in range(digits):
        pi += Decimal(1) / 16 ** k * (Decimal(4) / (8 * k + 1) - Decimal(2) / (8 * k + 4) - Decimal(1) / (8 * k + 5) - Decimal(1) / (8 * k + 6))
    return pi

# 求めたい円周率の桁数を入力
digits = int(input("計算したい円周率の桁数を入力してください: "))
pi = bbp_algorithm(digits)
print("円周率:", pi)

ニュートン・ラフソン（Newton-Paphson）法

参考サイト『Newton-Raphson法を用いた円周率近似』にあるとおり、次の漸化式を計算していきます。

$ \displaystyle L_{n+1} = \sqrt{2-\sqrt{4-(L_n)^2}}$

提示されたプログラム「ニュートン・ラフソン法」

ChatGPTから「ニュートン法」の名称と説明はありましたが、プログラムはアルゴリズム通りの処理ではなく、対応するプログラムは提示されませんでした。アルゴリズムの説明も適切ではありませんでした。

正多角形による方法（アルキメデスの方法）

円に内接・外接する正多角形を計算する方法6です。アルキメデスは、正 6 角形から出発して，辺の数が 2 倍にして，正 12 角形, 正 24 角形, 正 48 角形，正 96 角形，正 192 角形 について計算しました。

提示されたプログラム「正多角形による方法」

decimalライブラリを使ったプログラムが提示されいましたが、処理が正しく実装されていませんでした。

ヴィエトの公式

$\displaystyle \pi = 2 \cdot \frac{2}{\sqrt{2}} \cdot \frac{2}{\sqrt{2+\sqrt{2}}} \cdot \frac{2}{\sqrt{2+\sqrt{2+\sqrt{2}}}} \cdots$
の導出は正多角形によるものです。このヴィエトの公式7は ChatGPT では説明されませんでしたが、近しいプログラムが提案されましたので、取り上げました。

ラマヌジャン（Ramanujan）のアルゴリズム

$\displaystyle \frac{1}{\pi}= \frac{2\sqrt{2}}{99^{2}} \sum^{\infty}_{k=0} \frac{(4k)!(26390k+1103)}{(k!)^{4}396^{4k}}$

提示されたプログラム「ラマヌジャンのアルゴリズム」 
import math

def ramanujan(iterations):
    sum = 0
    for k in range(iterations):
        numerator = math.factorial(4 * k) * (1103 + 26390 * k)
        denominator = (math.factorial(k) ** 4) * (396 ** (4 * k))
        sum += numerator / denominator
    return 1 / ((2 * math.sqrt(2) / 9801) * sum)

# 収束の速さに応じて適切な反復回数を設定
iterations = int(input("計算したい円周率の桁数に応じた反復回数を入力してください: "))
pi = ramanujan(iterations)
print("円周率:", pi)

チュドノフスキー（Chudnovsky）の公式

$\displaystyle \frac{1}{\pi} =12\sum^{\infty}_{k=0} \frac{(-1)^{k}(6k)!(545140134k+13591409)}{(k!)^{3}(3k)!(-640320)^{3k+3/2}}$

提示されたプログラム「チュドノフスキーの公式」 
def chudnovsky(iterations):
    sum = 0
    for k in range(iterations):
        numerator = (-1) ** k * factorial(6 * k) * (545140134 * k + 13591409)
        denominator = factorial(3 * k) * factorial(k) ** 3 * 640320 ** (3 * k + 1.5)
        sum += numerator / denominator
    return 1 / (12 * sum)

def factorial(n):
    result = 1
    for i in range(1, n + 1):
        result *= i
    return result

# 収束の速さに応じて適切な反復回数を設定
iterations = int(input("計算したい円周率の桁数に応じた反復回数を入力してください: "))
pi = chudnovsky(iterations)
print("円周率:", pi)

算術幾何平均による方法

ガウス＝ルジャンドル（Gauss-Legendre）のアルゴリズム

$a_0 = 1, b_0 = \frac{1}{\sqrt{2}}$に設定して、以下の値を$n = 1, 2, 3, \dots$と順に計算します。
$\displaystyle a_{n + 1} = \frac{a_n + b_n}{2}$（$a_n$と$b_n$の相加平均）
$b_{n + 1} = \sqrt{a_n \cdot b_n}$（$a_n$と$b_n$の相乗平均）

$T = (a_0 - a_1)^2 + 2(a_1 - a_2)^2 + 4(a_2 - a_3)^2 + \cdots + 2^{n-1}(a_{n-1} - a_n)^2$

とすると、円周率$\displaystyle \pi = \frac{(a_n + b_n)^2}{1 - 4T}$と近似される。

提示されたプログラム「ガウス＝ルジャンドルのアルゴリズム」 
import math

def gauss_legendre(precision):
    a = 1.0
    b = 1.0 / math.sqrt(2)
    t = 1.0 / 4.0
    p = 1.0
    pi_old = 0

    for _ in range(precision):
        a_next = (a + b) / 2
        b = math.sqrt(a * b)
        t -= p * (a - a_next) ** 2
        a = a_next
        p *= 2
        pi = (a + b) ** 2 / (4 * t)
        
        # 収束判定
        if abs(pi - pi_old) < 10 ** -precision:
            break
        pi_old = pi

    return pi

# 任意の桁数で円周率を計算
num_digits = int(input("計算したい円周率の桁数を入力してください: "))
precision = int(math.log10(10 ** num_digits))
pi = gauss_legendre(precision)
print("円周率:", pi)

サラミン＝ブレント（Salamin–Brent）のアルゴリズム

$a_0 = 1, b_0 = \frac{1}{\sqrt{2}}, c_0^{2} = a_0^{2} - b_0^{2} $として、
$\displaystyle a_{n + 1} = \frac{a_n + b_n}{2}, b_{n + 1} = \sqrt{a_n b_n}, c_{n + 1} = \frac{1}{2}(a_n - b_n)$ とするとき、

$\displaystyle \pi = \lim_{n \to \infty} \frac{2 a_{n+1}^{2}}{1- \displaystyle \sum^{n}_{j=0} 2^{j} c_j^{2}}$

提示されたプログラム「サラミン＝ブレントのアルゴリズム」

ガウス＝ルジャンドルのアルゴリズムが提示されました。

ガンマ関数を計算して求める方法

$\displaystyle \Gamma(x)= \lim_{n \to \infty} \frac{n! \cdot n^{x} }{x(x+1)(x+2)\cdots(x+n)}$

を計算するアルゴリズムを用いて、$\displaystyle \Gamma(\frac{1}{2}) = \sqrt{\pi}$　を使って計算します。

提示されたプログラム「ガンマ関数を用いる方法」

$\displaystyle \Gamma(\frac{1}{4}) = \sqrt{2 \varpi \sqrt{2\pi}}$（$\varpi$はレムニスケート周率）を用いるプログラムが提示されました。

ゼータ関数を計算して求める方法

$\displaystyle \zeta(s)= \sum^{\infty}_{n=1} \frac{1}{n^{s}}$

を計算するアルゴリズムを用いて、$\displaystyle \zeta(2) = \frac{\pi^{2}}{6}$　を使って計算します。

提示されたプログラム「ゼータ関数を用いる方法」 
import math

def zeta_function(s, terms):
    result = 0.0
    for k in range(1, terms + 1):
        result += 1.0 / (k ** s)
    return result

def pi_from_zeta_function(terms):
    zeta_2 = zeta_function(2, terms)
    return math.sqrt(6 * zeta_2)

# 任意の桁数で円周率を計算
num_digits = int(input("計算したい円周率の桁数を入力してください: "))
terms = int(10 ** (num_digits / 2))
pi = pi_from_zeta_function(terms)
print("円周率:", pi)

ChatGPTの誤り

ChatGPTが誤っていた点で顕著だった内容を整理しました。

  • アルゴリズムとプログラムの相違：アルゴリズムの説明と提案されたプログラムが合っておらず、異なる手法で処理が書かれているケースがありました。
  • import mathを使う：そもそもmath.piの値を途中の計算で利用していたり、$\sqrt{2}$の値やルート根の計算をmathライブラリの関数で使った例が多く提示されました。
  • 用語の誤り：次のような間違いがありました。
    • ガウス＝アーコットの公式
      ”アーコット”は恐らくarccotのことだと推測します。arctanとarccotには関係式(B)が成り立つため、次式を提示しているものと思われます。
      $\displaystyle \begin{equation} \frac{\pi}{4}= 12 \arctan{\frac{1}{18}} + 8 \arctan{\frac{1}{57}} - 5 \arctan{\frac{1}{239}}\end{equation}$
    • チューリングの式
      Chudnovskyのアルゴリズムの一種として説明がされましたが、チューリングは円周率の計算分野には登場しなさそうです。
    • ボールウェインの公式
      "ボレルの加速法"という名称で説明がありましたが「ボールウェイン（Borwein）」であり、その説明には「リーマンのゼータ関数の無限和の値を使用」とあり、誤っているようです。
    • Brent-McMillanのアルゴリズム
      オイラー・マスケローニ定数$\displaystyle \gamma = \lim_{n \to \infty}(\sum^{n}_{k=1}\frac{1}{k} - \ln n) $を近似計算するプログラムが提示されました。

おわりに

ChatGPT が知っているアルゴリズムは以上になります。このほかに「モンテカルロ法」「連分数を用いる方法」等がありますが、回答には出てきませんでした。
本稿のキッカケになった100兆桁の円周率を計算したスクリプトは、GoogleがGitHubに公開しています。また、研究者の過去の動画がためになると思いました。最後にリンクを紹介して、本稿を締めることとします。

（●）（●） Happy Hacking!
/""　＿＿""＼

  1. Googleの研究者の記事『円周率の100兆桁目を求めて』 ↩

  2. 円周率.jp ↩

  3. 円周率公式のまとめ ↩

  4. 円周率を計算する ↩

  5. 高野喜久雄の公式
    $\displaystyle \begin{equation} \frac{\pi}{4}= 12 \arctan{\frac{1}{49}} + 32 \arctan{\frac{1}{57}} - 5 \arctan{\frac{1}{239}} + 12 \arctan{\frac{1}{110443}} \end{equation}$ ↩

  6. https://www.juen.ac.jp/math/nakagawa/openh20pi.pdf ↩

  7. https://www.kurims.kyoto-u.ac.jp/~ooura/pi04.pdf ↩

5
3
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
5
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?