fortran66のブログ

fortran について書きます。

OOP風エラトステネスの篩

OOP言語の常識に疎いので、コンストラクタ・デストラクタをどうするのが良いのかよく分かりません。デストラクタが二回呼ばれていますが、一回目はコンストラクタ内で仮生成したものが破棄されるところで出ているのかしらん?

デストラクタはF2003規格のfinal指定子によって(allocatableなら自動で解放されるのかもしれませんが・・)、コンストラクタはこの記事を参考にしました。“Object-Oriented Fortran: User-defined constructors”(これが良い方法なのかは不明です。)

またIntel Fortranでは、ユーザー定義I/Oがまだインプリメントされていないので、出力ルーチンは普通のサブルーチンで書きました。

実行結果

ソースプログラム

  module m_prime
    implicit none
    private
    public :: t_prime
    type :: t_prime
      integer, allocatable :: prm(:)
    contains
      procedure :: prn
      procedure :: prn_max
      final :: delete_prm
    end type t_prime
    
    interface t_prime ! consturctor alias
      module procedure init_prime
    end interface
    
  contains
  
    type (t_prime) function init_prime(np) ! constructor
      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 
      init_prime%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 function init_prime

    function eratos(n) result(res)
      integer, intent(in) :: n
      integer, allocatable :: res(:)
      logical(1) :: 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
    
    subroutine prn(this)
      class (t_prime), intent(in) :: this
      print '(10i8)', this%prm
      return
    end subroutine prn
    
    subroutine prn_max(this)
      class (t_prime), intent(in) :: this
      print *, this%prm( size(this%prm) )
      return
    end subroutine prn_max

    subroutine delete_prm(this) ! destructor
      type (t_prime), intent(in out) :: this
      if (allocated(this%prm)) deallocate(this%prm)
      print *, 'tiro finale! w'
      return
    end subroutine delete_prm
 
  end module m_prime
  
  program main
    use m_prime
    implicit none
    type (t_prime), allocatable :: prime
  ! Body of eratos
    prime = t_prime(10**2)
    call prime%prn()
    call prime%prn_max()
    deallocate(prime)
    stop
  end program main