C Check numbers for primality, using the precomputed list of
C primes that is generated by program PACKSV and is stored in
C the binary file PRIMP2G.BIN.
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C
C URL of this file:
C http://www.randomwalk.de/sequences/prim.for
C
C Version history:
C Sep 18 additional plausibility checks for user input,
C        improved error messages and output
C Sep 16 user input extende (e.g. suppress output of found primes
C        for negative input of upper search limit).
C Sep 3 initial version
      doubleprecision a, b, c
      integer n
      character line*80
      logical conly
C Initialize prime checking function, read stored packed table
      write (*,*) 'Start reading compressed table of primes.'
      call primin ( ipstat )
      if ( ipstat .ne. 0 ) then
        write (*,*) 'Error reading table, press Enter to continue'
        read (*,1000) line
        stop
      endif
C
C Start loop over user input lines
C
100   continue
C
C Ask user for input
      write (*,*)
     & 'Enter range to be checked,',
     & ' upper limit <0 to count only, 0 to stop'
C Target label for re-entering erroneous input
101   continue
      line = ' '
C read into character string
      read (*,1000) line
1000  format ( a )
C Find last non-blank character in line
      do 30 k=80,1,-1
C add terminating character for list-directed read
      if ( line(k:k) .ne. ' ' ) then
        line (k+1:k+1) = '/'
        goto 35
      endif
30    continue
35    continue
C Parse user input
      a = 1.0D0
      b = -1.25D0
      read (line,*,iostat=ios) a, b
      if ( ios .ne. 0 ) then
        write (*,*)
     & 'Numeric input required, please check and re-enter:'
        goto 101
      endif
C Check for termination
      if ( abs(a) .lt. 1.0D0 .and. abs(b) .lt. 2.0D0 ) stop
C If only one number was given, set tested interval to this number.
      if ( b .eq. -1.25D0 ) b = a
C Extract "count only flag" from second input number
      conly = b .lt. 0.0D0
C truncate negative input for lower limit
      if ( a .lt. 0.0D0 ) then
        a = 0.0D0
        write (*,*) 'Lower limit adjusted to 0.'
      endif
C Check range of user input
      b = abs(b)
      if ( a .gt. b ) then
        write (*,*)
     &  'Lower limit > upper limit. Nothing to check.'
        goto 100
      endif
      c = max(a,b)
C lower end of range
      if ( c .lt. 2.0D0 ) then
        write (*,*) 'Nothing to check.'
        goto 100
      endif
C deny testing ranges beyond 2>31
      if ( a .gt. 2.0D0**31-1.0D0 ) then
        write (*,*) 'Sorry, numbers>2^31-1 cannot be checked.'
        goto 100
      endif
C force upper limit
      if ( c .gt. 2.0D0**31-1.0D0 ) then
        write (*,*) 'Upper limit reduced to 2147483647'
        ib = 2147483647
      else
        ib = nint (c)
      endif
      ia = max ( 2, nint(a) )
C adjust search limit for even user input
      iaa = ia
      if ( mod(ia,2) .eq. 0 ) iaa = ia + 1
C
      if ( .not. conly ) write (*,*) 'Primes found:'
C catch special case: Lower boundary=2
      if ( ia .eq. 2 ) then
        n = 1
        if ( .not. conly ) write (*,*) 2
      else
        n = 0
      endif
C loop over odd numbers in requested range
      if ( conly ) then
        do 10 i = iaa, ib, 2
        if ( jsprim(i) ) n = n + 1
10      continue
C Output of result
        write (*,1002) ia, ib, n
1002  format ('Number of primes >=', I11, ' and <=', I11, ':', i9 )
      else
        do 11 i = iaa, ib, 2
        if ( jsprim(i) ) then
          write (*,*) i
          n = n + 1
        endif
11      continue
        if ( n .eq. 0 ) then
          write (*,*) 'None'
        else
          write (*,1002) ia, ib, n
        endif
      endif
C
C Skip back to ask user for next input
      goto 100
C
      end
C *********
      subroutine primin ( istat )
C Read and check packed prime table
C Author: Hugo Pfoertner http://www.pfoertner.org/
C
C Change history
C Sep 18 2004 added integrity checks
C Sep  3 2004 initial version
C
      parameter (ip=71582788)
      integer*1 pckdp(0:ip)
C Test btest(I1,I1) ->  integer*1 k
      common /packpr/ pckdp
C Preset status
      istat = 0
C Get packed prime table from unformatted file
      open ( unit=11, file='PRIMP2G.BIN', form='unformatted',
     &       status='old',iostat=ios )
      if ( ios .ne. 0 ) then
        write (*,1000)
1000  format ('Error reading PRIMP2G.BIN.',/, 
     & 'Must exist in the same directory, where prime.exe is called.',
     & /, 'Either copy PRIMP2G.BIN to present directory',/,
     & 'or create it by running packsv.exe.' )
        istat = 1
        goto 999
      endif
      read (11) pckdp
C Get signature
      ichksm = 0
      loc = 1
10    continue
      ichksm = ichksm + pckdp(loc)
      loc = loc + loc
      if ( loc .le. ip ) goto 10
      if ( ichksm .ne. 107 ) then
C Test   write (*,*) loc, ichksm
        write (*,*) 'Integrity check of PRIMP2G.BIN failed.'
        write (*,*) 'Re-create PRIMP2G.BIN by running packsv.exe'
        istat = 2
      endif
999   continue
      close (unit=11)
C End of subroutine PRIMIN
      end
C **********
      logical function jsprim ( z )
C Check for primality using compressed table of precomputed primes
C Author: Hugo Pfoertner hppt://www.pfoertner.org/
C
C Change history:
C Sep  3 2004 initial version
C
      implicit integer (a-z)
      parameter (ip=71582788,msmalp=31)
C
      integer k, map(0:29)
      integer*1 mk, pckdp(0:ip)
      common /packpr/ pckdp
      logical smlpri(0:msmalp)
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 smlpri /.false.,.false.,.true.,.true.,.false.,.true.,
     &                     .false.,.true.,.false.,.false.,.false.,
     &                     .true.,.false.,.true.,.false.,.false.,
     &                     .false.,.true.,.false.,.true.,.false.,
     &                     .false.,.false.,.true.,.false.,.false.,
     &                     .false.,.false.,.false.,.true.,.false.,
     &                     .true. /
C Start of test
C Special treatment for small numbers
      if ( z .le. msmalp ) then
C Look into table
        jsprim = smlpri(z)
      else
C determine location of relevant byte
        j = z / 30
C relative position in byte
        k = z - 30 * j
C Map to bits carrying information
        mk = map(k)
C standard locations for composites
        if ( mk .lt. 0 ) then
          jsprim = .false.
        else
C use pre-computed table values
          jsprim =  btest(pckdp(j),mk)
        endif
      endif
      return
C End of logical function JSPRIM
      end
