C Determine primes by sieving in a packed array that holds the
C primality information for the intervals (30*n,30*(n+1)) in one byte.
C Write the packed primality information for all numbers <2^31
C to a binary file with 71582788 bytes. The actual output file
C is slightly larger due to the header information generated by
C the Fortran unformatted write function.
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C
C Change history:
C Sep 16 2004 Ask user for progress indicator
C Sep 2 2004 Initial version
      implicit integer (a-z)
      parameter (ip=71582788,m=2147483647)
C
      dimension dif(0:7)
      integer*1 pckdp(0:ip), k, mk, map(0:29)
      character*1 c
C
      data map / -1,
     &   0, -1, -1, -1, -1, -1,  1, -1, -1, -1,
     &   2, -1,  3, -1, -1, -1,  4, -1,  5, -1,
     &  -1, -1,  6, -1, -1, -1, -1, -1,  7 /
      data dif / 1,  7, 11, 13, 17, 19, 23, 29 /
C
      mxsiev = 1 + nint ( sqrt ( float(m) ) )
      do 10 i = 0, ip
      pckdp(i) = -1
10    continue
      j = 0
      k = 0
C Start sieving with sieve width 7
      write (*,*) 'Start sieving. To watch progress, enter 1, else 0:'
      read (*,*) watch
      if ( watch .ne. 1 ) then
        write (*,*) 'Expected time to create table:'
        write (*,*) '270s / (GHz speed of processor)'
      endif
100   continue
      k = k + 1
      if ( k .gt. 7 ) then
        k = 0
        j = j + 1
      endif
      next = 30*j + dif(k)
      lim = m - next
      if ( btest(pckdp(j),k) ) then
C item is prime, use as sieve width
        if ( watch .eq. 1 ) write (*,1000) next
1000    format ( I6 )        
        d = next + next
200     continue
        jj = d / 30
        kk = d - 30 * jj
        mk = map(kk)
        if ( mk .ge. 0 ) THEN
           pckdp(jj) = ibclr ( pckdp(jj), mk )
        endif
        if ( d .le. lim ) then
          d = d + next
          goto 200
        endif
      endif
      if ( next .le. mxsiev ) goto 100
C Write to unformatted file
      open ( unit=11, file='PRIMP2G.BIN', form='unformatted',
     &       status='unknown' )
      write (*,*) 'Sieving finished, start writing to file.'
      write (11) pckdp
      close (unit=11)
C
      write (*,*) 'Primes written to binary file PRIMP2G.BIN'
999   continue
      write (*,*) 'To continue, press Enter.'
      read (*,1001) c
1001  format ( A ) 
C End of program PACKSV
      end
