GaussのLi(x)関数 - fortran66のブログ

Nまでの素数の個数の評価式として、ルジャンドルの $N\over \ln(N)$ やGaussの $Li(x)\equiv\int^N_2{1\over \ln(N)}$ などがあります。ここではガウスのLi(x)関数を数値積分で求めてみます。

定義式のままでは積分範囲が広すぎるので、 $y=\e x$ の変数変換を行い、さらに部分積分をしてみます。
$Li(x)=\int^{\ln N}_{\ln 2}{\e y\over y}dy=[\ln y \e y]^{\ln N}_{\ln 2}-\int^{\ln N}_{\ln 2}{\ln y \e y}dy$
これがベストなのか分かりませんが、なんとなく思いつきでｗ

N=10**3から10**11までの結果を求め、以下に示しました。分割区間数は最大2048までとってあります。おおむね収束しています。以下のプログラムでは、シンプソン法を呼び出すときには、分割区関数が偶数である必要があります。

台形積分を行うと、分割区間数を２倍にするごとに収束誤差が１／４づつに小さくなっているので、ちゃんと解析的な関数になっているようです。そういうわけで、リチャードソン補外をかけてやってよりよい収束値を得ることも出来ます。ここで調子に乗ってさらにもう何段かリチャードソン補外をかけてもいいのですが、それはやらずにシンプソン法でも数値積分を行うことにします。この時も、分割区間数を２倍にするごとに収束誤差が１／１６づつになっていることが確かめられたので、ここでもリチャードソン補外をかけてやります。リチャードソン補外をかけた値は、一番右側のカラムに表示しています。

得られた結果をネットの結果と比べると、0.5から1.0程度少なく出ていて謎です。定義どおりに計算していると思うのですが・・・・ま、今後の課題としておきますｗｗｗ

H22-2-18追記
理解しました。定義が違っていました。
一般的に、
$\mathrm{li}(x)\equiv\int^\infty_0{1\over \ln x}dx,$
$\mathrm{Li}(x)\equiv\int^\infty_2{1\over \ln x}dx=\mathrm{li}(x)-\mathrm{li}(2),$
となっていることが多いのですが、件のサイトでは、表記Li(x)で積分範囲の下限を0に取ったli(x)をあらわしているようです。li(2)=1.0451637801174927848445888891..なのでつじつまが合います。慣用的には二つの表記は混用されて使われているようです。(なおli(x)はx=1に特異点があるので主値積分で定義されていて普通に数値積分できません。）

結局、ここでの数値積分の結果も、上述のネットの結果も正しい値を与えておりました。表のすぐ上に定義式が書いてありした。（よく読んでなくてゴメン！ｗ）

Li(x) (the principal value of integral of 1/log u from u=0 to u=x)

ちなみに、多くの入門用書籍ではLi(x)とli(x)の定義が混同されていて、本文中の定義と表の値が約１ずれています。ご注意あれ。（今、手元にある本では４冊中２冊で間違っています。）
なお、RiemannのR(x)関数に対しては、どちらの定義で計算しても小数点以下２桁程度の誤差しか出ません。たしかメビウス関数の総和は０だったのでもっともな結果の気がします。

プログラムでは、公式をそのままFortranで表現するようにしたため無駄なテンポラリ・メモリーを暗黙に使っています。あまり伝統的なスタイルでありません。現代の計算機にとってはたいした問題じゃないんですが、らしくないですｗ

参照文献：『数値計算の常識』(1985), 伊理正夫, 藤野和建 .[主に'70年代の知見に基づく本ですが、面白い本です。]

実行結果

台形公式での結果

シンプソン法での結果

ソースプログラム

PROGRAM gauss
  IMPLICIT NONE
  INTEGER, PARAMETER :: kd = SELECTED_REAL_KIND(14), k = 1 ! 1:trapezoid,  2:simpson
  INTEGER :: i, j, n, i0 = 7, i1 = 11
  REAL(kd) :: x0, x1, d0, d1, y1, y0

  x0 = LOG(2.0_kd)
  d0 = 0.0_kd

  DO j = 3, 11
   x1 = LOG( 10.0_kd**j )
   y0 = 0.0_kd
   DO i = i0, i1
     n = 2**i
     y1 = f(x1) - f(x0) - trape(x0, x1, n) !simpson(x0, x1, n) !
     d1 = y1 - y0
     IF (i >= i0 + 2) &
       PRINT '(a, i2, i5, a, 2g15.6, g18.10)', 'x, n, Li(x): 10**', &
              j, n,':', y1, d0 / d1, (y1 - y0 / 4.0_kd**k) / (1.0_kd - 1.0_kd / 4.0_kd**k)
    d0 = d1
    y0 = y1
   END DO
   PRINT *
  END DO

  STOP

 CONTAINS

  REAL(kd) ELEMENTAL FUNCTION trape(x0, x1, n)
     REAL(kd), INTENT(IN) :: x0, x1
     INTEGER, INTENT(IN) :: n
     REAL(kd) :: h
     INTEGER :: i
     h = (x1 - x0) / n
     trape = h * ( 0.5_kd * ( f(x0) + f(x1) ) + SUM([( f(x0 + i * h), i = 1, n - 1 )])  )
     RETURN
  END FUNCTION trape

  REAL(kd) ELEMENTAL FUNCTION simpson(x0, x1, n)
     INTEGER, INTENT(IN) :: n
     REAL(kd), INTENT(IN) :: x0, x1
     REAL(kd) :: h, odd, even
     INTEGER :: i
     h = (x1 - x0) / n
     odd  = SUM([( f(x0 + i * h), i = 1, n    , 2  )])
     even = SUM([( f(x0 + i * h), i = 2, n - 2, 2  )])
     simpson = h / 3.0_kd * ( f(x0) + f(x1) + 4.0_kd * odd + 2.0_kd * even )
     RETURN
  END FUNCTION simpson

  REAL(kd) ELEMENTAL FUNCTION f(x)
     REAL(kd), INTENT(IN) :: x
     f = EXP(x) * LOG(x)
     RETURN
  END FUNCTION f

END PROGRAM gauss