エラトステネスの篩　ふたたびｗ

素数をｎ個もとめる。ガウスの素数定理に基づいて必要なエラトステネスの篩の大きさを推測する。求めたおおむねｎ個の素数の配列から最初のｎ個を取りだす。

オプションとして /assume:realloc_lhs　が必要。求める素数の数が大きい場合にはオプション /heap-arrays も必要となる。

実行結果

500までの素数の表

チェック
http://en.wikipedia.org/wiki/List_of_prime_numbers#The_first_500_prime_numbers

一応10**7の値でチェック　おｋ
http://primes.utm.edu/lists/small/millions/

原始プログラム

module m_prime
    implicit none
    integer, allocatable, save :: prm(:)
  contains
    subroutine init_prime(np)  
      integer, intent(in) :: np
      integer :: n
      ! np ~ n / log(n) ! approximate number of primes less than n : Gauss   
      ! approximate np-th prime n ~ np * log(n) ~ np * log( np * log(n) ) ...  
      n = np * log( np * log( real(np) ) ) + 3 
      prm = take(np, eratos(n))
      return
    contains
      function take(n, arr) result(res) 
        integer, intent(in) :: n, arr(:)
        integer, allocatable :: res(:)
        res = arr(1:n)
        return
      end function take
    end subroutine init_prime
  
    function eratos(n) result(res)
      integer, intent(in) :: n
      integer, allocatable :: res(:)
      logical :: table(n)
      integer :: i, j
      table = .true.
      table(1) = .false.
      do i = 2, sqrt( real(n) ) ! sieve of Eratosthenes
        if (table(i)) forall(j = 2 * i:n:i) table(j) = .false.
      end do
      res = pack([(i, i = 1, n)], table) ! filter [1,2,3..] by table:[F,T,T..]
      return
    end function eratos
    
  end module m_prime
  
  program main
    use m_prime
    implicit none
  ! Body of eratos
    call init_prime(500)
    print '(20i5)', prm
    stop
  end program main