fortran66のブログ

fortran について書きます。

添え字変数を使わない方向で素因数分解

DO LOOP のループ変数をできるだけ使わない、あるいは配列の添え字変数を使わない方向で、有限個の素数のテーブルにしたがって素因数分解し、テーブルにない素数は余りとして別変数に積を取ってためるプログラムを書いてみます。
試しなので、さほど一貫した方向性はありません。

実行結果

ソース・プログラム

module m_prime
  implicit none
  integer, parameter :: np = 10
  integer :: prime(np)
contains
  pure function mk_prime(n)
    integer, intent(in) :: n
    integer, allocatable :: mk_prime(:)
    integer :: k
    mk_prime = [2]
    k = 3
    do while(size(mk_prime) < n) 
      if ( all(mod(k, mk_prime) /= 0) ) mk_prime = [mk_prime, k]
      k = k + 2
    end do
    return
  end function mk_prime
  
  subroutine init_prime()
    prime = mk_prime(np)
    return
  end subroutine init_prime
end module m_prime

module m_pnum
  use m_prime
  implicit none
  type :: t_pnum
    integer :: numerator = 1, denominator = 1  
    integer :: p(np) = 0 
  end type t_pnum
contains
  pure elemental type(t_pnum) function to_pnum(n)
    integer, intent(in) :: n
    where (mod(n, prime) == 0) to_pnum%p = idivs(prime)
    to_pnum%numerator = n / product(prime**to_pnum%p, mod(n, prime) == 0)
    return
  contains
    pure elemental integer function idivs(k)
      integer, intent(in) :: k
      idivs = 1
      do
        if (mod(n, k**(idivs + 1)) /= 0) exit
        idivs = idivs + 1
      end do
      return
    end function idivs
  end function to_pnum
end module m_pnum

program main
  use m_pnum
  implicit none
  call init_prime()
  print *, pr_pnum([1:100]) !non-standard
  stop
contains
   pure elemental character(78) function pr_pnum(k)
     integer, intent(in) :: k
     write(pr_pnum, '(i6, "::", 2i6, ":", *(i4))') k, to_pnum(k) 
     return
   end function pr_pnum
end program main