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

はじめに

以前のエントリで、円周率 $\pi$ とネイピア数 $e$ を Binary Splitting 法という手法を使って、$\pi$ は1億桁まで、$e$ は1000万桁まで計算してみました。

本稿では、同じような手法で計算できる数『アペリー数』を取り上げ、Haskellのプログラムで計算してみましょう。

アペリーの定数

ゼータ関数を $\zeta$ とすると、$\zeta(3)$ で定義される数を『アペリーの定数（Apéry's constant）』といいます。フランスの数学者ロジェ・アペリーが、1978年に $\zeta(3)$ が無理数であることを示した『アペリーの定理』により命名されました。

リーマンゼータ関数は、$s$ を複素数、$n$ を自然数とするとき、

$$\begin{flalign}
\zeta(s) &:= \sum _{n=1} ^{\infty} \frac{1}{n^s} = 1 + \frac{1}{2^{s}} + \frac{1}{3^{s}} + \frac{1}{4^{s}} + \frac{1}{5^{s}} +\cdots &
\end{flalign}$$

で定義される関数 $\zeta$ のことをいいます。
この関数の $s=3$ の値がアペリーの定数です。

$$\begin{flalign}
\zeta(3) &= 1 + \frac{1}{2^{3}} + \frac{1}{3^{3}} + \frac{1}{4^{3}} + \frac{1}{5^{3}} +\cdots & \\
&= 1. 2020569031 5959428539 9738161511 4499907649 8629234049 \dots &
\end{flalign}$$

円周率.jpのBinary Splittingのページに$f(x)=2 \zeta(3) $のBinary Splitting Method（BS法）のための情報が、次のように紹介されています。

$f(x)$	$a(n)$	$b(n)$	$p(0)$	$q(0)$	$p(n)$	$q(n)$
$2\zeta(3)$	$(-1)^n (205n^2+250n+77) $	$1$	$1$	-	$n^5$	$32(2n+1)^5$

このための数式は Binary Splitting Recursion Libraryに記述があります。

$$\begin{flalign}
\zeta(3) &= \frac{1}{64}\sum _{k=0}^{\infty} \frac{(-1)^{k} (205k^{2}+250k+77)k! ^{10}}{(2k+1)! ^{5}} &
\end{flalign}$$

このときに、

$$\begin{flalign}
\zeta(3) &= \frac{T+77Q}{64Q} + O(1024^{-n}) &
\end{flalign}$$

となります。この数式を使って計算処理のプログラムを書くと次のようになります。

アペリーの定数 ζ(3) を求めるプログラム：CalcZ3.hs

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

main :: IO ()
main = print $ calcZ3 10_000_100

calcZ3 :: Integer -> Integer
calcZ3 nn = n10 * (t + 77*q) `div` (64*q)  where
  n10 = 10 ^ nn
  (_, _, q, t) = bisplit 0 nn (
      (\n->(if even n then 1 else -1)*(205*n*n+250*n+77))   -- an
    , (\n->1)     -- bn
    , (\n->n^5)   -- pn
    , (\n->32*(2*n+1)^5)   -- qn
    )

bisplit :: Integer -> Integer
  -> ((Integer -> Integer), (Integer -> Integer), (Integer -> Integer), (Integer -> Integer))
  -> (Integer, Integer, Integer, Integer)
bisplit nL nR fn
  | nR - nL == 1 =
    let (an, bn, pn, qn) = fn
    in (bn nR, pn nR, qn nR, (an nR) * (pn nR))
  | otherwise =
    let m = (nL + nR) `shiftR` 1
        (bL, pL, qL, tL) = bisplit nL m fn
        (bR, pR, qR, tR) = bisplit m nR fn
    in (bL * bR, pL * pR, qL * qR, bR * qR * tL + bL * pL * tR)

このプログラム（CalcZ3.hs）をghcでコンパイルし、処理時間を計測しながら実行してみます。なお、出力を標準出力にすると表示だけで時間がかかりますので、リダイレクトでファイルに出力します。1000万桁（プログラム上は誤差を考慮して＋100桁）までを計算するようにして結果をみてみましょう。

% ghc CalcZ3.hs
[1 of 2] Compiling Main ( CalcZ3.hs, CalcZ3.o )
[2 of 2] Linking CalcZ3
% time ./CalcZ3 > outputZ3
./CalcZ3 > outputZ3 219.87s user 27.61s system 84% cpu 4:53.21 total

この結果が正しいかをチェックするために、キリのよい桁数での値がZeta(3) - Apery's Constant ζ(3)に記載がありましたので、その値と比較してみます。コマンドで出力したアペリーの定数の値の1000万桁までの50個の数字は以下のとおり。

% tail -c151 outputZ3 | head -c50
65363268547751850215416264988233887058162791929459%

件のサイトに記載されている1000万桁までの直前の50個の数字は、

6536326854 7751850215 4162649882 3388705816 2791929459 : 10,000,000

ですので、出力された結果は正しそうです。
これ以外にも他の近似計算もありますが、参考サイトをご確認ください。¹

おわりに

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

本稿の環境

本稿のために使用した環境は以下となります。

macOS: Sonoma 14.6.1 (chip: Apple M1)
GHCup: 0.1.30.0
GHC: 9.6.6

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

https://en.wikipedia.org/wiki/Ap%C3%A9ry%27s_constant
https://archive.lib.msu.edu/crcmath/math/math/a/a274.htm
https://mathoverflow.net/questions/454929/ap%C3%A9rys-constant-zeta3-fastest-convergent-series ↩