エラトステネスのふるいで素数をしつこく数えます＋並列化
Fortran2003/08 では ビット操作関数が増えたので、これを用います。表は奇数のみ記録するとして、ある数 n に対して、配列の何番目の要素の何番目のビットかということを計算して、目的のビットを操作します。また、ある数 n が素数かどうかを返す論理型関数も定義します。
おおよそ 2*10^9 を超えると 32 ビットでは足りなくなってくるので、64ビット長でやっていますが、面倒なので行儀悪く手打ちで _8 を付けています。
実行結果
| 10^n | 10^n 以下の素数の個数 | 計算時間（秒） | 計算時間（並列）（秒） |
|---|---|---|---|
| 10^1 | 4 | ||
| 10^2 | 26 | ||
| 10^3 | 168 | ||
| 10^4 | 1,229 | ||
| 10^5 | 9,592 | ||
| 10^6 | 78,498 | ||
| 10^7 | 664,579 | ||
| 10^8 | 5,761,455 | 0.3 | 0.1 |
| 10^9 | 50,847,534 | 6 | 4.5 |
| 10^10 | 455,052,511 | 77 | 36 |
| 10^11 | 4,118,054,813 | 1643 | 396 |
記憶容量的には 10^12 まで行ける感じですが、10 時間以上かかりそうなのでやりませんでした。 無理っす、記憶容量足りません。計算違いしてましたｗ
実行時間は素数のテーブルの書き込みをする二重ループがほとんどを占めていて、そこは素数の数×テーブルの大きさ～lnN*N なので、この割合で実行時間が増えると思います。
ソース・プログラム
布団の中で考えたので適当です。
ビット操作関数としては ibitclr, ibtest, popcnt の三つを使いました。それぞれ、指定のビットを 0 にした整数値を返す、ビットが立っている時真を返す、立っているビットの数を返す関数です。
立っているビット数を数える sum(popcnt(itab)) が、10^1)~ の場合正しく働かず、32bit 演算をする感じなので、その時は DO LOOP を回して数えました。
program eratos_bit
implicit none
integer(8), parameter :: nmax = 10_8**1
integer(8) :: i, k, n, m, nb
integer(8), allocatable :: itab(:)
!
integer(8) :: ic0, icr0, ic1, icr1, ic2, icr2
call system_clock(ic0, icr0)
!
nb = 2_8 * bit_size(itab)
n = nmax / nb + 1
print *, 'array size=', n
allocate(itab(n))
itab = z'FFFFFFFFFFFFFFFF'!z'FFFFFFFF'
do i = nmax + 1, nb * n - 1, 2
call clr_bit(i)
end do
!
do i = 3, int(sqrt(real(nmax))), 2
if (is_prime(i)) then
do k = i**2, nmax, 2_8 * i
call clr_bit(k)
end do
end if
end do
!
call system_clock(ic1, icr1)
print *, 'eratos time', (ic1 - ic0) / real(icr1)
print *, sum(popcnt(itab))
! for 10^10~
!k = 0
!do i = 1, n
! k = k + popcnt(itab(i))
!end do
print *, 'no. of primes', k
call system_clock(ic2, icr2)
print *, 'bitcnt time', (ic2 - ic1) / real(icr2)
contains
subroutine clr_bit(n)
integer(8), intent(in) :: n
integer(8) :: ia, ib
ia = n / nb + 1_8
ib = mod(n, nb) / 2_8
itab(ia) = ibclr(itab(ia), ib)
end subroutine clr_bit
logical function is_prime(n)
integer(8), intent(in) :: n
integer(8) :: ia, ib
ia = n / nb + 1_8
ib = mod(n, nb) / 2_8
is_prime = btest(itab(ia), ib)
end function is_prime
end program eratos_bit
並列化
サブルーチンの部分を、ループ中に展開して、do concurrent と自動並列化で並列にしてみました。約４倍スピードアップ。
program eratos_bit
implicit none
integer(8), parameter :: nmax = 10_8**10
integer(8) :: i, k, n, m, nb
integer(8), allocatable :: itab(:)
!
nb = 2_8 * bit_size(itab)
n = nmax / nb + 1
print *, 'array size=', n
allocate(itab(n))
itab = -1_8 ! z'FFFFFFFFFFFFFFFF'!z'FFFFFFFF'
do i = nmax + 1, nb * n - 1, 2
block
integer(8) :: ia, ib
ia = i / nb + 1_8
ib = mod(i, nb) / 2_8
itab(ia) = ibclr(itab(ia), ib)
end block
end do
!
call stamp('begin of the sieve of Eratosthenes')
do i = 3, int(sqrt(real(nmax))), 2
if (is_prime(i)) then
do concurrent (k = i**2:nmax:2_8 * i)
block
integer(8) :: ia, ib
ia = k / nb + 1_8
ib = mod(k, nb) / 2_8
itab(ia) = ibclr(itab(ia), ib)
end block
end do
end if
end do
!
call stamp('end of the sieve of Eratosthenes')
!print *, sum(popcnt(itab))
! for 10^10 ~
k = 0
do i = 1, n
k = k + popcnt(itab(i))
end do
print *, 'no. of primes', k
call stamp('end of bit counting')
contains
logical function is_prime(n)
integer(8), intent(in) :: n
integer(8) :: ia, ib
ia = n / nb + 1_8
ib = mod(n, nb) / 2_8
is_prime = btest(itab(ia), ib)
end function is_prime
subroutine stamp(text)
character(*), intent(in) :: text
logical, save :: first = .true.
integer(8), save :: ic0, icr0
integer(8) :: ic1, icr1
if (first) then
call system_clock(ic0, icr0)
first = .false.
end if
call system_clock(ic1, icr1)
print '(f10.2, 2a)', (ic1 - ic0) / real(icr1), '(sec) : ', text
end subroutine stamp
end program eratos_bit