円周率を計算する（Haskell実装）

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

はじめに

少し前に『ChatGPTに円周率を計算するプログラムをお願いしてみた話』というエントリを投稿しました。ChatGPTが回答したプログラムは、mathパッケージやルートの計算（特に$\sqrt{2}$の値）をあてにしているプログラムが示されました。
本稿では、ライブラリの計算はあてにしないで円周率を計算するHaskellの実装を紹介します。最近はPCの性能もよく、次に挙げるように１億桁までは手元のPCで計算できるという記事も多くあります。

このitchynyさんのはてなブログでは10億桁に挑戦されている方も紹介されております。更には、Zennのtarohさんが100億桁に挑戦されている次のエントリもあります。

ここで紹介したエントリを投稿されている執筆者の方々の知識や熱意を紹介するだけで円周率を計算するための記事としては十分だと思いますが、改めてHaskellでの実装を試してみました。

円周率の計算

効率的な処理であるtanakhさんのGithubのHaskellプログラムを１つのプログラムにまとめて、紹介させて頂きます。

CalcPi.hs

CalcPi.hs

module Main (main) where
import Data.Bits (shiftR, shiftL)

main :: IO ()
main = print $ calcPI10 100_000_100

a, b, c, c3, c3div24 :: Integer
a = 13591409
b = 545140134
c = 640320
c3 = c * c * c
c3div24 = c3 `div` 24

calcPI10 :: Integer -> Integer
calcPI10 n = (pi2 * 10 ^ n) `shiftR` (fromIntegral n2) where
  n2 = floor (fromIntegral n * logBase 2 10 :: Double)
  pi2 = calcPI2 n2

-- Chudnovsky の公式をBinary Splitting 法で計算
calcPI2 :: Integer -> Integer
calcPI2 n = c3_2 * q `div` (12 * t + 12 * a * q) where
  c3_2 = c3 * (rsqrt2 n c3)
  (_, q, t) = bisplit 0 (n `div` 46)
  bisplit n1 n2
    | n2 - n1 == 1 =
      let an = (if even n2 then 1 else -1) * (a + b * n2)
          pn = (2 * n2 - 1) * (6 * n2 - 5) * (6 * n2 - 1)
          qn = n2 * n2 * n2 * c3div24
      in (pn, qn, an * pn)
    | otherwise =
      let m = (n1 + n2) `shiftR` 1
          (p1, q1, t1) = bisplit n1 m
          (p2, q2, t2) = bisplit m n2
      in (p1 * p2, q1 * q2, t1 * q2 + p1 * t2)

-- Newton法を用いた1/sqrt(n)の計算
rsqrt2 :: Integer -> Integer -> Integer
rsqrt2 d n
  | d < 8 = floor $ approx * 2 ^ d
  | otherwise = iter 8 $ (floor $ approx * 2 ^ initd) `shiftL` (di - initd)
  where
    di = fromInteger d
    initd = min 300 di
    approx = 1.0 / sqrt (fromIntegral n :: Double)   
    iter !p !x
      | p >= d = x
      | otherwise = iter (p*2) ((3 * x - n * x `mult` x `mult` x) `shiftR` 1) where
        mult x y = (x * y) `shiftR` di

このプログラムをghcでコンパイルし、処理時間を計測しながら実行してみます。出力を標準出力にすると表示だけで時間がかかりますので、リダイレクトでファイル出力にします。

% ghc CalcPi.hs
[1 of 2] Compiling Main ( CalcPi.hs, CalcPi.o )
[2 of 2] Linking CalcPi
% time ./CalcPi > outputPi
./CalcPi > outputPi 258.83s user 9.05s system 99% cpu 4:30.47 total

4分30秒程で1億桁まで計算できたようです。
前出の『円周率を1億桁計算しました！ ― その試行錯誤の詳しい経緯と結果 ー』のエントリで『一億桁を｢きちんと｣超えておく』の欄で検証されている値を参考にして、1億桁を超えた近辺の結果を確認してみましょう。円周率小数点以下1億1番目の数は2, 次は1, そして5... となり、10文字ずつにすると次の値なのだそうです。

$\begin{aligned}
& 2150588095\ 7832796348\ 7309513528\ 4911033417\ 9757201258\\
& 8340621369\ 0542295838\ 789{\color{red}4}607142\ 4855972210\ 0848156605
\end{aligned}$

今回紹介したプログラムでは、1億100桁まで出力していますので、最後の方（以下の出力）をこの値と比較してみます。すると、1億73桁までは一致しますが、赤色のところで数値が一致していません。そして、それ以降は一致していないことが分かります。

% tail -c 101 outputPi
21505880957832796348730951352849110334179757201258
83406213690542295838789582519940916909387067982357

このことから、指定した桁数よりも少し少ない桁までしか正確に計算できていないことになります。つまり、このプログラムでは目的の桁数より少し多めの桁まで計算しなければならないことが分かります。上の結果から、1億100桁までを指定して計算すれば、1億桁までは正確に計算できていると考えてよいでしょう。

Chudnovsky の公式

それではプログラムの処理の概略を解説します。まず、元となる$\pi$の計算式は、Chudnovsky の公式です。

Chudnovsky の公式

$\begin{eqnarray} \frac{1}{\pi} &=& \frac{12}{\sqrt{C^3}}\sum_{n=0}^{\infty} \frac{(-1)^n(6n)!(A+Bn)}{(3n)!(n!)^3C^{3n}}\\ \\ A &=& 13591409 = 13\cdot1045493\\ B &=& 545140134 = 2\cdot3^2\cdot7\cdot11\cdot19\cdot127\cdot163\\ C &=& 640320 = 2^6\cdot3\cdot5\cdot23\cdot29 \end{eqnarray}$

これを円周率.jpのChudnovsky の公式の式変形にあるように式変形します。

$\displaystyle K_n = \frac{(-1)^n(6n)!}{(3n)!(n!)^3C^{3n}}$とおきかえると、$\displaystyle \frac{1}{\pi} = \frac{12}{\sqrt{C^3}}\sum_{n=0}^{\infty} K_n(A+Bn)$となります。そして次の関係式が成り立ちます。
$$\begin{eqnarray}
\frac{K_{n}}{K_{n-1}} &=&
\frac{-1\cdot(6n)(6n-1)(6n-2)(6n-3)(6n-4)(6n-5)}
{(3n)(3n-1)(3n-2)n^3 C^3}\\
&=&\frac{-(6n-1)(2n-1)(6n-5)}{n^3 (C^3/24) }\tag{A}
\end{eqnarray}$$

この式を用いてBS法で計算します。BS法については円周率.jpのBinary Splittingに解説がありますが、下記に説明を載せておきます。

Binary Splitting Method（BS法）

BS 法は、求めようとする級数が
$$\begin{aligned}
S = \sum_{n=0}^{\infty} \frac{a(n)}{b(n)} \frac{p(0)\cdots p(n)}{q(0)\cdots q(n)}
\end{aligned}$$の形で表わされているとします。このとき、この部分級数
$$\begin{aligned}
S = \sum_{n_1 \lt n\lt n_2}
\frac{a(n)}{b(n)} \frac{p(n_1)\cdots p(n)}{q(n_1)\cdots q(n)}
\end{aligned}$$を直接計算するのではなく、$P=p(n_1)\cdots p(n_2-1)$、$Q=q(n_1)\cdots q(n_2-1)$、$B=b(n_1)\cdots b(n_2-1)$を考えて、$T=BQS$を求めてから$S=T/BQ$を導く方法です。

$m=\lfloor(n_1+n_2)/2\rfloor$として、$n_1\leqq n\lt m$に属する項を$P_l, Q_l, B_l, T_l$、$m\leqq n\lt n_2$に属する項を$P_r, Q_r, B_r, T_r$とするとき、
$$B = B_lB_r,\ P = P_lP_r,\ Q = Q_lQ_r,\ T = B_rQ_rT_l + B_lP_lT_r$$で再帰的に計算することができます。

プログラムCalcPi.hsでは、式（A）を使い、次の関数を当てはめてBS法で計算しています。
an = (if even n2 then 1 else -1) * (a + b * n2)
pn = (2 * n2 - 1) * (6 * n2 - 5) * (6 * n2 - 1)
qn = n2 * n2 * n2 * c3div24
（プログラム中には明記されませんが、bn = 1としています）

平方根の計算

次に、Chudnovsky の公式では、$C = 640320$としたときの$\sqrt{C^{3}}$を用います。そこで平方根を高精度に計算する必要があります。平方根の計算ロジックは、円周率.jpのNewton法にも詳しくありますが、ニュートン法を用いて$1/\sqrt{A}$を計算し、その値を用いて$\sqrt{A}=A \times 1/\sqrt{A}$で計算します。

プログラムCalcPi.hsでは、関数rsqrt2 :: Integer -> Integer -> Integerが逆平方根を計算する関数で、$\sqrt{C^{3}}$を計算している箇所はc3_2 = c3 * (rsqrt2 n c3)です。

おわりに

ほかにもシフト演算を活用したりして計算の工夫がされていますが、計算に大きく関連しているのは「$\pi$を求める公式」「BS法」「逆平方根の計算とニュートン法」をおさえておくとよいでしょう。

お手元の環境で検証する際に、動作が異なるときには参考になるかもしれませんので、最後に、本稿を記載するために検証したHaskell環境を記しておきます。

本稿の環境

本稿のために使用した環境は以下となります。
macOS: Sonoma 14.4 (chip: Apple M1)
GHCup: 0.1.22.0
GHC: 9.6.4

本稿は、tanakhさんのプログラムを掲載して一部を解説したエントリです。また、前述した記事もかなり参考にしております。tanakhさんはもちろんのこと、各記事の執筆者の皆さまに感謝いたします。

また、ご一読いただきまして有り難うございます。
（●）（●） Happy Hacking!
/""　＿＿""＼