fortran66のブログ

fortran について書きます。

エラトステネスの篩

確かに抽象度の高い新しい言語使用を使うとソースコードは簡潔になるけれど、見えないところで動的にメモリーを使用するので、スタックがあふれたり、最適化が効かなかったりで、DO LOOPのさわやかさには及びません。


出発ベクトルを演算子で何段も処理してゆくのを関数でやろうとして、苦労した時のことを思い出すと、合成演算子としてのポイントフリーや、中途のパイプのテンポラリ排除の為の遅延評価、帰り値にエラーを含ませるMaybe型が、大体どういう動機で湧いてきたのか透けて見える気がします。 

実行結果

ソースコード

無駄な[1..nmax]を生成するので、すぐスタックがあふれます。

PROGRAM primes
   IMPLICIT NONE
   INTEGER, PARAMETER   :: nmax = 1000
   INTEGER, ALLOCATABLE :: iprime(:)

   CALL eratos()
   PRINT *, iprime
   PRINT *, 'nmax =', nmax, 'number of primes =', SIZE( iprime )
   
   STOP

CONTAINS
   
   SUBROUTINE eratos()
      LOGICAL, ALLOCATABLE :: table(:) 
      INTEGER :: i, j
! make table    
      ALLOCATE( table(nmax) )
      table(1)  = .FALSE.
      table(2:) = .TRUE.
      DO i = 2, nmax
       IF ( table(i) ) FORALL (j = 2 * i:nmax:i) table(j) = .FALSE.
      END DO
! pick up prime numbers 
      ALLOCATE( iprime(COUNT( table )) )
      iprime = PACK( [(i, i = 1, nmax)], table ) 
      DEALLOCATE(table)
    
      RETURN
   END SUBROUTINE eratos 
  
END PROGRAM primes

その他

10**10までの素数を求めた場合。

n番 ルジャンドルの式 k/ln(k) n番目の素数 (k)

ソース

PROGRAM test
   IMPLICIT NONE
   INTEGER(8), PARAMETER :: nmax = 10_8**3 !2_8**33
   INTEGER(8), ALLOCATABLE :: iprime(:)

   CALL eratos()
   PRINT *, SIZE( iprime )
   PRINT *, iprime
   
   STOP

CONTAINS
   
   SUBROUTINE eratos()
      LOGICAL(1), ALLOCATABLE :: table(:) 
      INTEGER(8) :: i, j
! make table    
      ALLOCATE( table(nmax) )
      table = .TRUE.
      table(1) = .FALSE.
      DO i = 1, nmax
       IF ( table(i) ) FORALL (j = 2 * i:nmax:i) table(j) = .FALSE.
      END DO
! pick up prime numbers 
      ALLOCATE( iprime(COUNT( table )) )
       
      i = 0
      DO j = 1, nmax
        IF (table(j)) THEN 
          i = i + 1
          iprime(i) = j
        END IF
      END DO   
      DEALLOCATE(table)
    
      RETURN
   END SUBROUTINE eratos 
  
END PROGRAM test