0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@kyamaz

Binary Splitting法とその例 - ζ(5)の数値計算 - （Haskell実装）

0
Last updated at Posted at 2026-10-31

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

はじめに

以前に投稿した @kyamaz :frog:のエントリで、円周率 $ \pi $ の1億桁までの計算や ネイピア数 $ e $ の1000万桁までの計算を紹介しました。その計算では、Binary Splitting Method という手法をご紹介しました。具体的な数値計算のプログラムもHaskellの実装例として取り上げました。

また更に、リーマンゼータ関数 $\zeta(s) $ で $ s = 3 $ の値、アペリーの定数 と呼ばれる $ \zeta(3) $ の計算も紹介しました。

本稿は、Binary Splitting 法を用いた計算例の第４弾になります。

ゼータ関数の奇数の値

さて突然ですが、リーマンゼータ関数で最も興味を惹かれているのは「リーマン予想」だと思います。ただ、リーマンゼータ関数はそれ以外にも未解決な問題が知られています。引数 $s$ が整数のケースに限ってその偶奇で分かっていることと分かっていないことの特徴をみてみましょう。（$ n \in \mathbb{N}$のとき）

  • 偶数の値 $\zeta(2n)$ は、ベルヌーイ数 $B_{2n}$ を使って次のように表されることが分かっており、完全に解明されています。つまり、$\pi$ を用いて簡潔に表せるということになります。

    $\displaystyle \zeta(2n)=(-1)^{n+1}\frac{B_{2n}(2\pi)^{2n}}{2(2n)!} $

  • 一方で、奇数の値 $\zeta(2n+1)$ に関しては、ほとんど分かっていません。

    • $\zeta(3)$ が無理数であることは、アペリーにより示されていますが、$\zeta(3)$ が超越数であるかどうかは知られていません。
    • $\zeta(5), \zeta(7), \dots $ が有理数なのか無理数なのかすらもわかっていません。
    • $\zeta(2n)$ のように $\pi$ や有理数の組み合わせで表す公式は現状では発見されていません。多くの高速収束級数（Ramanujan 型、Koecher 型）が知られていますが、解析的に単純な公式は知られていません。（閉形式で簡単に書けません。）

などなど。未解決な問題も多いです。
そのなかでも、研究されている数学者の方もいて、それなりに進展しているようです。例えば、

  • 「$\zeta(5), \zeta(7),\zeta(9), \zeta(11)$ のうち少なくとも1つは無理数」と示されています。（Zudilin, 2001）2
  • 「奇数 $\zeta(2n+1)$ の無理数が無限に存在する」ことが示されています。（Rivoal, 2000）3

「無理数になる奇数ゼータ値の個数」の下限は、その後も Fischler–Sprang–Zudilin（2019）4 や Lai（2025）5 によって改良が続いています。また、2024年には Calegari–Dimitrov–Tang が Apéry の方法を大きく広げ、$L(2,\chi_{-3})$ などの無理性を示しました67。ただし、$\zeta(5)$ 単独の無理性は今も未解決です。

Binary Splitting 法

前出のエントリでも紹介しましたが、Binary Splitting 法（BS法）は、級数で表された計算を、分割して合算していく手法です。区間 $ [a,b) $ の部分和を有理数 $P(a,b)/Q(a,b)$ として再帰的に作りあげます。左右の部分和をそれぞれ $P_L/Q_L, P_R/Q_R$ としたときに合成ルールとして、$ P = P_L \cdot Q_R + P_R \cdot Q_L, Q = Q_L \cdot Q_R $ を使って深さ $ O(\log{N})$の木までまとめていく手法です。割り算を最小にし、乗算（大きな整数）を階層的にまとめることで、全体の演算コストを大幅に下げられる利点があります。なお、中間での計算は、通常、整数／有理数で扱われるため丸め誤差も生じません。

BS 法は、特に超幾何級数（hypergeometric series）型で展開された計算に有効です。典型的には次のような形：

$$\begin{flalign}
S = \sum_{n=0}^{\infty} \frac{a(n)}{b(n)}
\end{flalign}$$

ここで $a(n), b(n)$ は多項式あるいは階乗を含む整数列です。多くの場合 $b(n)$ は階乗の積や冪乗などで「積として書ける」形のときに有効な計算手法となります。

さて、少し数式を使って説明しましょう。

通常の計算では、$n$ の項ごとに分数を足していきますが、Binary Splitting では「区間和を分子・分母にまとめて持つ」という考え方をとります。
区間 $ [m,n) $ の部分和を
$$\begin{flalign}
S(m,n) = \sum_{k=m}^{n-1} \frac{p(k)}{q(k)}
\end{flalign}$$ とすると、これを
$$\begin{flalign}
S(m,n) = \frac{P(m,n)}{Q(m,n)}
\end{flalign}$$ という形にすることを狙います。ここで、$P(m,n), Q(m,n)$ は共に整数とします。

再帰的合成則

区間 $ [m,n) $ を真ん中で分割して $[m,l),[l,n) $ にすると、明らかに
$$\begin{flalign}
S(m,n) = S(m,l) + S(l,n)
\end{flalign}$$ となり、これを分数で書くと、
$$\begin{flalign}
\frac{P(m,n)}{Q(m,n)} = \frac{P(m,l)}{Q(m,l)} + \frac{P(l,n)}{Q(l,n)}
\end{flalign}$$ 通分して、左辺と右辺を比べると、
$$\begin{flalign}
P(m,n) &= P(m,l)Q(l,n) + P(l,n)Q(m,l) \\
Q(m,n) &= Q(m,l)Q(l,n)
\end{flalign}$$ これが Binary Splitting の核心となる合成則 です。

つまり、「和を分母の積で通分してまとめる操作」を、区間分割で構造化したものにすぎません。この利点は、区間分割することで、大きな積を一度に作るのではなく、部分積をバランスよく組み合わせることで、計算の「桁数の成長」に対して効率的になります。

さて、前出のエントリで紹介した BS 法はこの考え方を発展させて、求めようとする級数が
$$\begin{flalign}
S = \sum_{n=0}^{\infty} \frac{a(n)}{b(n)} \frac{p(0) \cdots p(n)}{q(0) \cdots q(n)}
\end{flalign}$$ の形で表されているとします。このとき、区間 $n_1 \lt n \leqq n_2$ の部分級数を、区間の先頭から積を取り直した形
$$\begin{flalign}
S(n_1,n_2) = \sum_{n_1 \lt n \leqq n_2}
\frac{a(n)}{b(n)} \frac{p(n_1+1) \cdots p(n)}{q(n_1+1)\cdots q(n)}
\end{flalign}$$ で考えます。これを直接計算するのではなく、

$P=p(n_1+1) \cdots p(n_2),\ Q=q(n_1+1) \cdots q(n_2),\ B=b(n_1+1) \cdots b(n_2)$ を考えて、$T=BQ,S(n_1,n_2)$ を求めてから $S(n_1,n_2)=T/BQ$ を導く方法として一般化しました。

項が1つだけ（$n_2 = n_1 + 1$、$n = n_2$）のときは、$B = b(n),\ P = p(n),\ Q = q(n),\ T = a(n)p(n)$ です。
$m=\lfloor(n_1+n_2)/2\rfloor$ として、$n_1\lt n\leqq m$に属する項を$P_l, Q_l, B_l, T_l$、$m\lt n\leqq 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$$で再帰的に計算することができます。

Haskellの実装具体例は以下のようになります。

import Data.Bits (shiftR)

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)

ζ(5)の数値計算

まずは比較のために、定義そのままの級数

$\displaystyle \zeta (5) = \sum _{n=1}^{\infty} \frac{1}{n^5} $

を使って計算してみましょう。各項を有理数として扱い、$ t_{n} = 1/n^5 = p_{n} / q_{n} $ とみて、$ p_{n} = 1, q_{n} = n^5 $ とおきます。前の節の最初に説明した、単純な合成則（$P, Q$ だけを持つ形）をそのまま使います。

{-# LANGUAGE BangPatterns #-}
-- Binary splitting for zeta(5) using exact rationals (Integer arithmetic)

import Data.Ratio
import Numeric (showFFloat)

-- single term: p_n / q_n  with p_n = 1, q_n = n^5
termPQ :: Integer -> (Integer, Integer)
termPQ n = (1, n ^ 5)

-- bisplit over [l .. r)  (r exclusive)
-- returns (P, Q) representing sum_{n=l}^{r-1} p_n/q_n = P / Q
bisplit :: Integer -> Integer -> (Integer, Integer)
bisplit l r
  | l >= r    = (0,1)                     -- empty sum
  | l + 1 == r = termPQ l                 -- single term
  | otherwise =
      let m = (l + r) `div` 2
          (pL,qL) = bisplit l m
          (pR,qR) = bisplit m r
          -- combine: P = pL*qR + pR*qL, Q = qL*qR
          !p = pL * qR + pR * qL
          !q = qL * qR
      in (p, q)

-- wrapper: compute zeta(5) partial sum with N terms (n=1..N)
zeta5PartialRational :: Integer -> Rational
zeta5PartialRational n
  | n <= 0    = 0 % 1
  | otherwise = let (p,q) = bisplit 1 (n+1) in p % q

-- tail bound: for zeta(s) with s>1, integral test gives
-- R_N = sum_{n=N+1}^\infty 1/n^s <= integral_N^\infty x^{-s} dx = 1/((s-1) N^{s-1})
-- For s=5: R_N <= 1/(4 * N^4)
tailBoundForZeta5 :: Integer -> Double
tailBoundForZeta5 n = 1.0 / (4.0 * fromIntegral (n ^ (4 :: Int)))

-- convenience: print Rational to decimal with given digits (uses Double, limited precision)
rationalToDecimal :: Rational -> Int -> String
rationalToDecimal r digs =
  let d = fromRational r :: Double
  in showFFloat (Just digs) d ""

-- example main: compute partial sum N and show estimate
main :: IO ()
main = do
  let digits = 12           -- digits to display (Double-limited)
      -- choose N by tail bound to get ~digits of accuracy:
      -- want 1/(4 N^4) <= 10^{-digits}  => N >= (10^{digits} / 4)^{1/4}
      targetDigits = 12
      nGuess = ceiling ((10 ** fromIntegral targetDigits / 4) ** (1/4 :: Double)) :: Integer
      n = max 100 nGuess   -- at least some minimal N
  putStrLn $ "Using N = " ++ show n ++ " terms (estimated tail <= " ++ show (tailBoundForZeta5 n) ++ ")"
  let r = zeta5PartialRational n
  putStrLn $ "Rational sum = " ++ show (numerator r) ++ " / " ++ show (denominator r)
  putStrLn $ "Approx (Double, " ++ show digits ++ " digits) = " ++ rationalToDecimal r digits
  putStrLn $ "Tail bound <= " ++ show (tailBoundForZeta5 n)

実行例

Using N = 708 terms (estimated tail <= 9.949631079319851e-13)
Rational sum = 6087481254040475413754318371110047254112213895970965897653... / ...（分子・分母とも約1500桁なので省略）
Approx (Double, 12 digits) = 1.036927755142
Tail bound <= 9.949631079319851e-13

最も素直な $\displaystyle \sum _{n=1}^{\infty} \frac{1}{n^5} $ を使いました。これは実装が単純な反面、収束は遅めです。708 項も足して、ようやく12桁程度（誤差は $1/(4N^4) \approx 10^{-12}$ 以下）です。また、各項 $1/n^5$ は前の節で述べた「積として書ける」形にはなっていないので、BS 法の利点も十分には生かせていません。より速い方法（少ない項数で同じ精度）として Ramanujan / Koecher 型の高速級数を使う方がよいです。

Ramanujan / Koecher 型

Borwein & Bradley が引用している Koecher の式の形：

$$
\zeta(5) = 2 \sum _{k=1}^{\infty} \frac{(-1)^{k+1}}{k^5 \pmatrix{2k \\ k}} - \frac{5}{2} \sum _{k \ge 1} \frac{(-1)^{k+1} H _{k-1} ^{(2)}}{k^3 \pmatrix{2k \\ k}}
$$

ここで $\displaystyle H _{k-1} ^{(2)} := \sum _{j=1} ^{k-1} \frac{1}{j^2}$ （２次調和数）です。

この式は (A) 中央二項係数 $\displaystyle \pmatrix{ 2k \\ k } $ が分母に現れること、(B) 第２項目は $ H _{k-1} ^{(2)} $（累積和）を含む点が特徴です。この式を使う利点は、分母の $\displaystyle \pmatrix{ 2k \\ k } $ がおよそ $4^k$ の速さで大きくなるため、項が $4^{-k}$ 程度で小さくなり、非常に早く収束する点です（1項あたり約0.6桁ずつ精度が上がります）。さらに、隣り合う項の二項係数の比が $\displaystyle \pmatrix{2k \\ k} \Big/ \pmatrix{2k-2 \\ k-1} = \frac{2(2k-1)}{k}$ と簡単な有理式になるので、形が hypergeometric / Apéry 型の構造を持っており Binary-Splitting に馴染む点です。

Binary-Splitting に落とし込む方針（要点）

Koecher 式の RHS を
$$
S _{1} := \sum _{k \ge 1} \frac{(-1)^{k+1}}{k^5 \pmatrix{2k \\ k}},\quad
S _{2} := \sum _{k \ge 1} \frac{(-1)^{k+1} H _{k-1} ^{(2)}}{k^3 \pmatrix{2k \\ k}}
$$

として $\displaystyle \zeta(5) = 2 S_{1} - \frac{5}{2} S_{2} $ とします。$S_1, S_2$ それぞれを、単純な合成則の BS 法で計算します。

import Data.Ratio

-- Koecher-type zeta(5) with Binary Splitting
-- We also add a decimal expansion printer for Rational numbers.

-- central binomial coefficient C(2k,k)
centralBinom :: Integer -> Integer
centralBinom k
  | k <= 0    = 1
  | otherwise = product [(k+1)..(2*k)] `div` product [1..k]

-- precompute H2[k] = sum_{j=1}^k 1/j^2 as Rationals for k=0..N
precomputeH2 :: Integer -> [Rational]
precomputeH2 n = scanl (\acc j -> acc + 1 % (j*j)) 0 [1..n]

-- term for S1: (-1)^{k+1} / (k^5 C(2k,k))
termS1 :: Integer -> (Integer, Integer)
termS1 k =
  let sign = if odd k then 1 else -1
      num  = sign
      den  = k^5 * centralBinom k
  in (num, den)

-- term for S2: (-1)^{k+1} H_{k-1}^{(2)} / (k^3 C(2k,k))
termS2 :: [Rational] -> Integer -> (Integer, Integer)
termS2 h2 k =
  let sign = if odd k then 1 else -1
      c    = centralBinom k
      h    = if k <= 1 then 0 else h2 !! fromIntegral (k-1)
      hn   = numerator h
      hd   = denominator h
      num  = sign * hn
      den  = hd * k^3 * c
  in (num, den)

-- Binary Splitting
bisplitSeq :: (Integer -> (Integer,Integer)) -> Integer -> Integer -> (Integer,Integer)
bisplitSeq term l r
  | l >= r     = (0,1)
  | l + 1 == r = term l
  | otherwise  =
      let m = (l + r) `div` 2
          (pL,qL) = bisplitSeq term l m
          (pR,qR) = bisplitSeq term m r
          p = pL * qR + pR * qL
          q = qL * qR
      in (p,q)

-- partial sums
computeS1S2 :: Integer -> (Rational, Rational)
computeS1S2 n
  | n <= 0    = (0,0)
  | otherwise =
      let h2 = precomputeH2 n
          (p1,q1) = bisplitSeq termS1 1 (n+1)
          (p2,q2) = bisplitSeq (termS2 h2) 1 (n+1)
      in (p1 % q1, p2 % q2)

zeta5Approx :: Integer -> Rational
zeta5Approx n =
  let (s1,s2) = computeS1S2 n
  in 2 * s1 - (5 % 2) * s2

---------------------------------------------------------
-- Decimal expansion of a Rational, arbitrary precision
---------------------------------------------------------
showDecimal :: Int -> Rational -> String
showDecimal digits r =
  let (q,rn) = quotRem (numerator r) (denominator r)
      intPart = show q
      fracDigits = go rn (denominator r) digits
  in intPart ++ "." ++ concat fracDigits
  where
    go :: Integer -> Integer -> Int -> [String]
    go _ _ 0 = []
    go n d k =
      let (q,rn) = quotRem (n*10) d
      in show q : go rn d (k-1)

---------------------------------------------------------
-- Example run
---------------------------------------------------------
main :: IO ()
main = do
  putStrLn "Approximate zeta(5) using Koecher series (Binary Splitting)"
  let n = 200   -- number of terms
      r = zeta5Approx n
  putStrLn $ "With N=" ++ show n ++ " terms:"
  putStrLn $ "zeta(5) ≈ " ++ showDecimal 50 r

実行例（N=200, 50 桁表示）

Approximate zeta(5) using Koecher series (Binary Splitting)
With N=200 terms:
zeta(5) ≈ 1.03692775514336992633136548645703416805708091950191

N=200 項で、表示した50桁はすべて正しい値になっています。素朴な級数では708項で12桁程度だったのと比べると、収束の速さの違いがよく分かります。

おわりに

いかがでしたでしょうか？Binary Splitting 法（BS法）の特徴を生かした計算方法を $\zeta(5)$ を題材にしてみてみました。気になる定数をこの計算手法で試してみるのも面白いかもしれません。
最後に、本稿を記載するために検証したHaskell環境を記しておきます。お手元の環境で検証する際に、動作が異なるときには参考になるかもしれません。

本稿の環境

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

macOS: Tahoe 26.6.2 (chip: Apple M1)
GHCup: v0.2.6.2
GHC: 9.12.2

ご一読いただきまして有り難うございます。

（●）（●） Happy Hacking!
/""　＿＿""＼

  1. @kyamaz は、オープンソース・コミュニティ『OpenQL』プロジェクト8を通じて、皆さんと共に量子情報・量子コンピューティングの分野で挑戦しております。引き続きどうぞ宜しくお願い致します。 ↩

  2. W. Zudilin, "One of the numbers ζ(5), ζ(7), ζ(9), ζ(11) is irrational", Russian Mathematical Surveys 56:4 (2001), 774–776 — https://www.mathnet.ru/eng/rm427 ↩

  3. T. Rivoal, "La fonction zêta de Riemann prend une infinité de valeurs irrationnelles aux entiers impairs", C. R. Acad. Sci. Paris, Sér. I 331 (2000), no. 4, 267–270. arXiv:math/0008051 — https://arxiv.org/abs/math/0008051 ↩

  4. S. Fischler, J. Sprang, W. Zudilin, "Many odd zeta values are irrational", Compositio Mathematica 155 (2019), 938–952. arXiv:1803.08905 — https://arxiv.org/abs/1803.08905 ↩

  5. L. Lai, "A note on the number of irrational odd zeta values, II", arXiv:2501.05321 (2025). $\zeta(3), \zeta(5), \dots, \zeta(s-1)$ のうち少なくとも $1.284\sqrt{s/\log s}$ 個が無理数であることを示したもの — https://arxiv.org/abs/2501.05321 ↩

  6. F. Calegari, V. Dimitrov, Y. Tang, "The linear independence of $1$, $\zeta(2)$, and $L(2,\chi_{-3})$", arXiv:2408.15403 (2024) — https://arxiv.org/abs/2408.15403 ↩

  7. E. Klarreich, "Rational or Not? This Basic Math Question Took Decades To Answer", Quanta Magazine, 2025-01-08 — https://www.quantamagazine.org/rational-or-not-this-basic-math-question-took-decades-to-answer-20250108/ ↩

  8. OpenQLプロジェクトは、量子コンピューターを扱うためのライブラリを開発するためのオープンソースプロジェクトです。量子情報、量子コンピューターに興味のある人たちが集うコミュニティを運営しております。詳しくはconnpassのサイトをご覧ください。 ↩

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?