fortran66のブログ

fortran について書きます。

素数の個数の近似関数

Riemannによる素数の個数の式。
\mathrm{R}(x)\equiv \sum^\infty_{k=1}(\mu(k)/k)\mathrm{Li}(k^{1/k}).

ここで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