夏休みFortran祭り 素数を10^11まで数える

#エラトステネスのふるいで素数をしつこく数えます＋並列化

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 を回して数えました。

```fortran
    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 と自動並列化で並列にしてみました。約４倍スピードアップ。

```fortran
    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
```