円周率の計算：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}^{\mathrm{\infty }}\frac{(-1{)}^m(6m)!(A+Bm)}{(3m)!(m!{)}^3{C}^{3m}}}}}$$

$$\begin{array}{r}\{\begin{array}{l}A=13591409\\ B=545140134\\ C=640320\end{array}\end{array}$$

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 -