Riemannによる素数の個数の式。
ここでLi(x)はGaussによる素数の個数の近似式です。(でも積分の下限が2なのか0なのかよく分からんw ただ計算結果からは、どちらでも大差無し)またμ(k)はメビウス関数です。
ソースプログラム
あまりまじめに考えなかった。
PROGRAM Riemann IMPLICIT NONE INTEGER, PARAMETER :: kd = SELECTED_REAL_KIND(14) INTEGER :: k DO k = 1, 15 PRINT '(a, i2, 2F20.2)', ' 10**', k, Riemann_R(10.0_kd**k) END DO STOP CONTAINS REAL(kd) FUNCTION riemann_R(x) REAL(kd), INTENT(IN) :: x INTEGER, PARAMETER :: moebius(60) = & [ 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1, 0, 1, 1, -1, 0, 0, 1, 0, 0, -1, -1, & -1, 0, 1, 1, 1, 0, -1, 1, 1, 0, -1, -1, -1, 0, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, 1, 0, 1, 1, -1, 0 ] INTEGER :: k REAL(kd) :: r(60) r = [( moebius(k) * gauss_Li(10, x**(1.0_kd / REAL(k, kd)), trapezoid), k = 1, 60 )] riemann_r = SUM( [( r(k) / REAL(k, kd), k = 1, 60 )] ) RETURN END FUNCTION riemann_R REAL(kd) FUNCTION gauss_Li(n, x, func) INTEGER, INTENT(IN) :: n REAL(kd), INTENT(IN) :: x INTERFACE REAL(kd) FUNCTION func(x0, x1, k) IMPORT REAL(kd), INTENT(IN) :: x0, x1 INTEGER, INTENT(IN) :: k END FUNCTION func END INTERFACE REAL(kd) :: x0, x1 REAL(kd), ALLOCATABLE :: y(:) INTEGER :: i, j x0 = LOG( 2.0_kd) x1 = LOG( x ) y = f(x1) - f(x0) - [( func(x0, x1, 2**i), i = 1, n )] DO j = 1, n - 1 y = [( richardson(y(i), y(i + 1), j), i = 1, SIZE(y) - 1 )] END DO gauss_Li = y(1) RETURN END FUNCTION gauss_Li REAL(kd) ELEMENTAL FUNCTION richardson(x0, x1, k) REAL(kd), INTENT(IN) :: x0, x1 INTEGER, INTENT(IN) :: k richardson = (x1 - x0 / 2.0_kd**k) / (1.0_kd - 1.0_kd / 2.0_kd**k) RETURN END FUNCTION richardson REAL(kd) ELEMENTAL FUNCTION trapezoid(x0, x1, n) REAL(kd), INTENT(IN) :: x0, x1 INTEGER, INTENT(IN) :: n REAL(kd) :: h INTEGER :: i h = (x1 - x0) / n trapezoid = h * ( 0.5_kd * ( f(x0) + f(x1) ) + SUM([( f(x0 + i * h), i = 1, n - 1 )]) ) RETURN END FUNCTION trapezoid 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 Riemann