円周率の計算：Chudnovskyの公式

• 1 概要
• 2 公式
• 3 Chudnovskyの公式による計算
• 4 参考文献

1 概要

チュドノフスキーアルゴリズム は、ラマヌジャン系の円周率を求める公式で \(\pi\) の桁を計算するための高速な方法です。1988年にChudnovsky兄弟によって発表され、2021年8月14日に62.8兆桁の世界記録計算に使用されました。

2 公式

\[ \displaystyle \frac{1}{\pi} = \displaystyle \frac{12}{ \sqrt{C^3}} \displaystyle \sum_{m=0}^{\infty} \frac{(-1)^m(6m)!(A+Bm)}{(3m)!(m!)^3C^{3m}} \]

\[ \begin{eqnarray} \left\{ \begin{array}{l} A=13591409 \\ B=545140134 \\ C=640320 \end{array} \right. \end{eqnarray} \]

3 Chudnovskyの公式による計算

library(Rmpfr)

chudo <- function(n, prec) {
    A <- mpfr(13591409, prec)
    B <- mpfr(545140134, prec)
    C <- mpfr(640320, prec)
    S <- mpfr(0, prec)

    for (m in 0:n) {
        mm <- mpfr(m, prec)
        f6m <- factorial(6*mm)
        f3m <- factorial(3*mm)
        fm3 <- factorial(mm)^3
        S <- S + (-1)^m * f6m * (A + B * m) / (f3m * fm3 * C^(3 * m))
    }

    return (C^(3/2) / (12 * S))
}

digit <- 150
prec <- 5 + ceiling(digit / log10(2))
ConstPi <- Const('pi', prec)

m <- 0:9
p <- sapplyMpfr(m, chudo, prec)

tibble(m, digit=-log10(as.numeric(abs(ConstPi-p)))) %>%
    ggplot() +
    aes(x=m, y=digit) +
    geom_col() +
    geom_text(aes(label=sprintf('%4.1f', digit)), color='white', vjust=1.1)

\(m\) が1進む毎に精度が約14桁上がっている。

4 参考文献

• Ramanujan系公式 - 円周率.jp -
• Chudnovsky の公式を用いた円周率の計算用メモ - Qiita -
• Chudnovsky algorithm - Wikipedia -
• 円周率 - Wikipedia -