ゼータ関数の数値計算（Haskell実装）

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

はじめに

今どきは、Vibe Coding と呼ばれていますが、AIによるコーディング支援も活用すれば、単機能なプログラムを作成するエントリは流行らないだろうと思っていますが、ゼータ関数を数値計算するためのHaskellの実装という少々マニアック（？）なプログラムは、即座に動作するプログラムには行きつかなかったため、記事にしてみます。

リーマン予想を検証したい？

リーマン予想といれば、リーマンゼータ関数の全てのゼロ点は $ z = \frac{1}{2} $ に乗っているという予想ですが、ゼータ関数の数値計算で確認したくなるのが世の常です。（←限られた興味の持ち主の気持ち）
先人たちも様々な方法で調べていたり、プログラムの実装例があります。

これに習って、Haskellで実装してみようというのが本稿の主旨です。

Haskellで精度を気にした計算をするには、MPFRを使うのが良さそうです。MPFR (Multiple-Precision Floating-Point Reliably) は、高精度な浮動小数点演算を行うためのC言語ライブラリです。GNU MPライブラリ (GMP) を基盤としており、標準の倍精度浮動小数点数よりも桁数の多い任意精度で計算できることが特徴です。正確な丸め処理を重視し、科学技術計算や正確な数値計算に利用されます。
Haskellでは、hmpfrというパッケージがあり、内部でMPFRのC言語ライブラリを使えます。プログラムを書く前に、このパッケージのインストールをしておきましょう。

事前準備

必要な C ライブラリ（libgmp-dev および libmpfr-dev）を、OS のパッケージマネージャーでインストールします。

Debian/Ubuntu

sudo apt-get install libgmp-dev libmpfr-dev

macOS

brew install gmp mpfr

stack/cabal でパッケージをインストール（パスの指定はmacOSの例）

stack build --extra-include-dirs=/opt/homebrew/include --extra-lib-dirs=/opt/homebrew/lib hmpfr

cabal install --lib hmpfr \
  --extra-include-dirs=/opt/homebrew/include \
  --extra-lib-dirs=/opt/homebrew/lib

ゼータ関数の数値計算の実装

Data.Number.MPFRにはPreludeと同じ関数が定義されているため注意が必要です。

Zeta.hsのソースコード

Zeta.hs

{-# LANGUAGE FlexibleContexts #-}

module Zeta (C, zeta, showC, fromDoubleC) where

import Prelude hiding (sqrt, log, exp, sin, cos, atan2, div)
import Data.Number.MPFR
  ( MPFR
  , Precision
  , fromDouble
  , RoundMode(Near)
  , sqrt, log, exp, sin, cos, atan2
  , add, sub, mul, div
  , toString
  )
import Data.Complex

-- 型エイリアス: 複素数 (MPFR)
type C = Complex MPFR

-- 整数から MPFRへの変換
fromIntegerPrec :: Precision -> Integer -> MPFR
fromIntegerPrec prec n = fromDouble Near prec (fromIntegral n)

-- Double から複素数（MPFR）への変換
fromDoubleC :: Precision -> Double -> Double -> C
fromDoubleC prec re im = fromDouble Near prec re :+ fromDouble Near prec im

-- 複素数べき乗 (x ^ y)
expC :: Precision -> C -> C -> C
expC prec (a :+ b) (c :+ d) =
  let aa = mul Near prec a a
      bb = mul Near prec b b
      r = sqrt Near prec (add Near prec aa bb)
      theta = atan2 Near prec b a
      ln_r = log Near prec r
      c_ln_r = mul Near prec c ln_r
      d_theta = mul Near prec d theta
      exp_arg = sub Near prec c_ln_r d_theta
      exp_real = exp Near prec exp_arg
      d_ln_r = mul Near prec d ln_r
      c_theta = mul Near prec c theta
      exp_imag = add Near prec d_ln_r c_theta
      cos_imag = cos Near prec exp_imag
      sin_imag = sin Near prec exp_imag
      realPart = mul Near prec exp_real cos_imag
      imagPart = mul Near prec exp_real sin_imag
  in realPart :+ imagPart

-- k^-s
powC :: Precision -> Int -> C -> C
powC prec k s =
  let kMPFR = fromIntegerPrec prec (toInteger k)
      zeroMPFR = fromIntegerPrec prec 0
  in expC prec (kMPFR :+ zeroMPFR) s

-- 複素数の除算
divC :: Precision -> C -> C -> C
divC prec (a :+ b) (c :+ d) =
  let cc = mul Near prec c c
      dd = mul Near prec d d
      denom = add Near prec cc dd
      ac = mul Near prec a c
      bd = mul Near prec b d
      bc = mul Near prec b c
      ad = mul Near prec a d
      realNum = add Near prec ac bd
      imagNum = sub Near prec bc ad
      realPart = div Near prec realNum denom
      imagPart = div Near prec imagNum denom
  in realPart :+ imagPart

-- 複素数の加算
addC :: Precision -> C -> C -> C
addC prec (a :+ b) (c :+ d) = add Near prec a c :+ add Near prec b d

-- 複素数の減算
subC :: Precision -> C -> C -> C
subC prec (a :+ b) (c :+ d) = sub Near prec a c :+ sub Near prec b d

-- Dirichlet η 部分和
etaPartial :: Precision -> C -> Int -> C
etaPartial prec s maxN = go 1 (fromIntegerPrec prec 0 :+ fromIntegerPrec prec 0)
  where
    go k acc
      | k > maxN = acc
      | otherwise =
          let sign = if odd k then 1 else -1
              term = divC prec (fromIntegerPrec prec (toInteger sign) :+ fromIntegerPrec prec 0) (powC prec k s)
          in go (k + 1) (addC prec acc term)

-- 複素数の絶対値
magnitudeC :: Precision -> C -> MPFR
magnitudeC prec (a :+ b) =
  let aa = mul Near prec a a
      bb = mul Near prec b b
  in sqrt Near prec (add Near prec aa bb)

-- ζ(s) = η(s) / (1 - 2^(1-s))
-- 精度 prec, 部分和 maxN
zeta :: Precision -> C -> Int -> Maybe C
zeta prec s maxN
  | nearOne   = Nothing
  | otherwise = Just $ divC prec (etaPartial prec s maxN) (subC prec one (expC prec two (subC prec one s)))
  where
    one = fromIntegerPrec prec 1 :+ fromIntegerPrec prec 0
    two = fromIntegerPrec prec 2 :+ fromIntegerPrec prec 0
    s_minus_one = subC prec s one
    mag = magnitudeC prec s_minus_one
    epsilon = fromDouble Near prec 1e-30
    nearOne = let magStr = toString 10 mag
              in read magStr < 1e-30

-- 複素数を文字列表示 (小数 digits 桁)
showC :: C -> Int -> String
showC (a :+ b) digits =
  let digitsW = fromIntegral digits :: Word
      showMPFR m = toString digitsW m
  in showMPFR a ++ " + " ++ showMPFR b ++ "i"

ビルドは以下のようにパッケージを指定します。

Zeta.hsのビルド

ghc -package hmpfr Zeta.hs

プログラムから呼び出すには、次の実装ファイルをビルドして、実行します。

Main.hs

import Zeta

main :: IO ()
main = do
  let prec = 200     -- 計算精度（ビット, 200 ≈ 60桁）
      maxN = 200000  -- η 部分和の項数
      s = fromDoubleC prec 0.5 14.1347251417346937904572519835625

  let Just v = zeta prec s maxN
  putStrLn (showC v 40)

また、GHCiからは以下のように呼び出せます。

ghci>   let prec = 200     -- 計算精度（ビット, 200 ≈ 60桁）
ghci>       maxN = 200000  -- η 部分和の項数
ghci>       s = fromDoubleC prec 0.5 14.1347251417346937904572519835625
ghci>
ghci>   let Just v = zeta prec s maxN
ghci>   putStrLn (showC v 40)

どちらの場合も出力は次のようになります。

出力結果

4.184938781034103006207381999880932313e-4 + 2.161444339696390041031618744137440338e-4i

なお、下記サイトに、ゼロ点のリストがあります。その最初のゼロ点を目掛けて検証してみましたが、結果はそこそこの精度のようです。

おわりに

Vibe Coding では自分の環境に合った実装が提案されませんでしたので、記事にしてみましたが、このような記事の価値はほぼ無くなってしまっているかもしれませんね。

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

本稿の環境

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

macOS: Sequoia 15.6.1 (chip: Apple M1)
GHCup: 0.1.50.2
GHC: 9.6.7

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