Problem 012

i番目の三角数は i(i+1)/2 素因数分解は i, i+1, 2 についてやって足し引きする。約数の数は素因数に1足して積。

ソース・プログラム

    program PEuler012
      implicit none
      integer :: i, k1, k2, n
      integer, allocatable :: ipfac(:), ipfac0(:), ipfac1(:)
      i = 1
      ipfac0 = [integer::]
      do
        ipfac1 = ipfactor(i + 1)
        k1 = size(ipfac0)
        k2 = size(ipfac1)
        if (k1 > k2) then 
          ipfac = ipfac0
          ipfac(:k2) = ipfac(:k2) + ipfac1
        else
          ipfac = ipfac1
          ipfac(:k1) = ipfac(:k1) + ipfac0
        end if
        ipfac(1) = ipfac(1) - 1
        n = product( ipfac + 1, ipfac /= 0 )
        if (n > 500) exit
        i = i + 1
        ipfac0 = ipfac1
      end do
      print *, i, ':', i * (i + 1) / 2, ':', n
      stop
    contains
      function ipfactor(n)
        integer, value :: n
        integer, allocatable :: ipfactor(:)
        integer, allocatable :: iptab(:), jp(:)
        iptab = ieratos(n)
        ipfactor = spread(0, 1, ncopies = size(iptab))
        do
          jp = mod(n, iptab)
          where (jp == 0) ipfactor = ipfactor + 1
          n = n / product(iptab, jp == 0)
          if (n == 1) exit
        end do
        return
      end function ipfactor
    
      function ieratos(n) ! thieve of Eratostenes
        integer, intent(in) :: n
        integer, allocatable :: ieratos(:)
        logical :: tab(n)
        integer :: i
        tab = .true.
        tab(1) = .false.
        do i = 2, int( sqrt(real(n)) ) 
          if ( tab(i) ) tab(i * i:n:i) = .false.
        end do
        ieratos = pack( [ (i, i = 1, n) ], mask = tab )
        return
      end function ieratos      
    end program PEuler012

出力結果

12375 : 76576500 : 576
続行するには何かキーを押してください . . .