Nまでの素数の個数の評価式として、ルジャンドルのやGaussの
などがあります。ここではガウスのLi(x)関数を数値積分で求めてみます。
定義式のままでは積分範囲が広すぎるので、の変数変換を行い、さらに部分積分をしてみます。
これがベストなのか分かりませんが、なんとなく思いつきでw
N=10**3から10**11までの結果を求め、以下に示しました。分割区間数は最大2048までとってあります。おおむね収束しています。以下のプログラムでは、シンプソン法を呼び出すときには、分割区関数が偶数である必要があります。
台形積分を行うと、分割区間数を2倍にするごとに収束誤差が1/4づつに小さくなっているので、ちゃんと解析的な関数になっているようです。そういうわけで、リチャードソン補外をかけてやってよりよい収束値を得ることも出来ます。ここで調子に乗ってさらにもう何段かリチャードソン補外をかけてもいいのですが、それはやらずにシンプソン法でも数値積分を行うことにします。この時も、分割区間数を2倍にするごとに収束誤差が1/16づつになっていることが確かめられたので、ここでもリチャードソン補外をかけてやります。リチャードソン補外をかけた値は、一番右側のカラムに表示しています。
得られた結果をネットの結果と比べると、0.5から1.0程度少なく出ていて謎です。定義どおりに計算していると思うのですが・・・・ま、今後の課題としておきますwww
H22-2-18追記
理解しました。定義が違っていました。
一般的に、
となっていることが多いのですが、件のサイトでは、表記Li(x)で積分範囲の下限を0に取ったli(x)をあらわしているようです。li(2)=1.0451637801174927848445888891..なのでつじつまが合います。慣用的には二つの表記は混用されて使われているようです。(なおli(x)はx=1に特異点があるので主値積分で定義されていて普通に数値積分できません。)
結局、ここでの数値積分の結果も、上述のネットの結果も正しい値を与えておりました。表のすぐ上に定義式が書いてありした。(よく読んでなくてゴメン!w)
Li(x) (the principal value of integral of 1/log u from u=0 to u=x)
ちなみに、多くの入門用書籍ではLi(x)とli(x)の定義が混同されていて、本文中の定義と表の値が約1ずれています。ご注意あれ。(今、手元にある本では4冊中2冊で間違っています。)
なお、RiemannのR(x)関数に対しては、どちらの定義で計算しても小数点以下2桁程度の誤差しか出ません。たしかメビウス関数の総和は0だったのでもっともな結果の気がします。
プログラムでは、公式をそのままFortranで表現するようにしたため無駄なテンポラリ・メモリーを暗黙に使っています。あまり伝統的なスタイルでありません。現代の計算機にとってはたいした問題じゃないんですが、らしくないですw
参照文献:『数値計算の常識』(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