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

はじめに

前に『円周率を計算する（Haskell実装）』のエントリを投稿いたしました。そこで、円周率$\pi$を求める計算をされている猛者の方々の記事を紹介しましたが、そのなかでtarohさんが下に紹介する記事で『$e$ (自然対数の底) を100億桁求める』ことに挑戦されております。

円周率が計算できるようになると、他の無理数（特に超越数）も計算したくなるものです。そこで本稿では、tarohさんにならって、$e$ (自然対数の底) を計算するプログラムを取り上げてます。Haskellで実装した円周率を計算するプログラムが手元にありますので、ネイピア数もHaskellで実装してみましょう。

ネイピア数

WikiPediaのネイピア数によれば、ネイピア数は、古くは1618年ジョン・ネイピアによって発表された対数の研究の付録に収録されていた近似値の表があり、自然対数の底がネイピア数と呼ばれているのはここから来ており、記号として$e$が用いられるのは、数学者レオンハルト・オイラーにちなんでつけられているようです。
定義は、

$\displaystyle
e = \lim_{n \to +\infty} \Big( 1 + \frac{1}{n} \Big) ^{n}
$

です。また、$f(x)=e^{x}$のマクローリン展開から次式も成り立ちます。

$\displaystyle
e = \sum_{n=0}^{\infty} \frac{1}{n!}
$

ネイピア数を計算する

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

$f(x)$	$a(n)$	$b(n)$	$p(0)$	$q(0)$	$p(n)$	$q(n)$
$\exp(x)$	1	1	1	1	u	nv

※表中の$x$は$x=u/v$

ここで、$u=1, v=1$として実装すれば、ネイピア数$e$を計算する処理が書けそうです。

Binary Splitting Method

Binary Splitting Method（BS法）の処理を少し汎用的にして、Haskellで実装してみましょう。第三引数のfnは4つのラムダ式の組(an,bn,pn,qn)をとるようにしてあります。

Binary Splitting法

import Data.Bits (shiftR, shiftL)
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)

e を計算する

前述の表から、$a(n)=1, b(n)=1, p(n)=1, q(n)=n$ を使えば、目的の値を計算する処理になります。

ネイピア数 e を計算する

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

calcE :: Integer -> Integer
calcE nn = n10 + (t * n10) `div` q  where
  n10 = 10 ^ nn
  (_, _, q, t) = bisplit 0 nn (
      (\n->1)   -- an
    , (\n->1)   -- bn
    , (\n->1)   -- pn
    , (\n->n)   -- qn
    )

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

% ghc CalcE.hs
[1 of 2] Compiling Main ( CalcE.hs, CalcE.o )
[2 of 2] Linking CalcE
% time ./CalcE > outputE
./CalcE > outputE 25.84s user 6.22s system 90% cpu 35.490 total

さて、結果が気になるところです。そこで1000万桁のデータをダウンロードできるサイトがありましたので、そこから$e$の値をダウンロードして、計算した結果と比較してみましょう。

tail -c 100 e10_000_000.txt # 正解の値
45724092905469358878328181949178345076018520370348
94317635589923616102851120105544429298561396705376
% tail -c 201 outputE # 計算した値
45724092905469358878328181949178345076018520370348
94317635589923616102851120105544429298561396705376
16484089932201833912492779274270196740747544061929
77002991551759476095731736714496572432599234884647

双方を比較すると、赤字のところまで一致しており、正しい値が計算できているようです。

おわりに

もう一桁大きな桁数まで計算しようとして1億桁に挑戦したところ、わたしの環境ではメモリ不足になりましたので、本稿では1000万桁までの計算とさせて頂きました。

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

本稿の環境

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

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