Sieve of Eratosthenes
julia と fortran でおおむね同じなプログラムで比較。10**9 以下の素数をエラトステネスの篩で求める。
Julia 1.5.1
julia さん速くて Fortran -O3 と揺らぎの範囲で大差ない。
@time pr = prim(10^9) 33.018580 seconds (5 allocations: 14.901 GiB, 0.30% gc time)
function prim(n::Int64) tab = Vector(1:n) tab[1] = 0 for i = 2:floor(Int64, sqrt(n)) if tab[i] != 0 for j = i*i:i:n tab[j] = 0 end end end filter(!iszero, tab) end
50847534-element Array{Int64,1}: 2 3 5 7 11 13 17 19 23 29 31 37 41 ⋮ 999999667 999999677 999999733 999999739 999999751 999999757 999999761 999999797 999999883 999999893 999999929 999999937
Intel Fortran 19.1.2
/heap-arrays100000 として大きい配列を stack から heap に移さないと stack overflow する。
int64 50847534 33.68400 int32 50847534 27.19900
今の時代の計算機はメモリー移動が律速段階だから 32bit にすると速い。もしくは Fortran の default 長が 32bit のためかもしれない。64 bit 整数対応は後付け。
module prim_m use, intrinsic :: iso_fortran_env implicit none contains function prime(n) result(ires) integer, intent(in) :: n integer(int64), allocatable :: ires(:) integer(int64) :: i, itab(n) forall (i = 1:n) itab(i) = i itab(1) = 0 do i = 2, int(sqrt(real(n))) if (itab(i) /= 0) itab(i*i::i) = 0 end do ires = pack(itab, itab /= 0) end function prime end module prim_m program prim use prim_m implicit none integer(int64), allocatable :: ip(:) integer(int64) :: ic0, ic1, irate call system_clock(ic0) ip = prime(10**9) call system_clock(ic1, irate) print *, size(ip), (ic1 - ic0) / real(irate) end program prim
数値計算のためのFortran90/95プログラミング入門(第2版)
- 作者:牛島 省
- 発売日: 2020/01/28
- メディア: 単行本(ソフトカバー)