■

よくあるFortranの例題として、エラトステネスのふるいを用いた素数計算があります。
これを用いて素数の数を求め、ガウスによる素数定理を調べることにします。

ここでエラトステネスのふるいとは、ヘレニズム時代のアレキサンドリアの自然学者エラトステネスによるといわれる素数を求めるアルゴリズムで、まず２以上ある数Nまでの正の整数のリストを作ります。次にリストの最小の数2を素数として記録し2の倍数をリストから消します。同様にリスト上の最小の数3を素数として記録しその倍数をリストから消します。以下5,7・・・について同様の操作を行います。倍数を消す操作は√Nまでの数についてやれば十分です。

ガウスによる素数定理とはｎ以下の素数の数がで近似できるというものです。

一億以下の素数の数で見ると大体6%程度の誤差で、近似が成り立っていることが分かります。

PROGRAM eratos
IMPLICIT NONE
INTEGER, PARAMETER :: np = 8, nmax = 10**np
LOGICAL(1) :: table(nmax)
INTEGER :: i, j, n, ix
REAL :: x
! Sieve of Eratosthenes
n = INT(SQRT(REAL(nmax))) 
table = .FALSE.
table(1) = .TRUE.
DO i = 2, n 
 IF (table(i)) CYCLE 
 DO j = i**2, nmax, i
  table(j) = .TRUE.
 END DO
END DO
! Primenumber Theorem
PRINT *, 'the number of prime numbers below 10**i'
DO i = 1, np
 ix = COUNT(.NOT. table(1:10**i))
 x  = 10.0**i / LOG(10.0**i)
 PRINT *, i, ix, x, ix / x
END DO
STOP
END PROGRAM eratos 

プログラムでは、素数の表は１バイトのLOGICAL変数を用いて、非素数を『真』、素数を『偽』として表しました。Fortran90らしい命令としては、配列の『真』要素を数えるCOUNT()関数を用いたことでしょうか。

fortran66.hatenablog.com