C----------------------------------------------------------------------
C SIEVE7.FOR
C
C Simple implementation of Sieve of Eratosthenes search for primes.
C
C 18/04/2020
C----------------------------------------------------------------------
C Version with direct loop limits to avoid overflow
C----------------------------------------------------------------------
      SUBROUTINE FINDP5(nkfac, kfac, kmax, kprimes, nprimes, kkmax,
     +           iFAIL)
      IMPLICIT NONE
      INTEGER*4 nkfac               ! vector length
      INTEGER*4 kfac(1:nkfac)       ! vector of already-known primes
      INTEGER*4 kmax                ! vector length
      INTEGER*4 kprimes(1:kmax)
      INTEGER*4 nprimes             ! always less than kmax
      INTEGER*4 kkmax
      INTEGER*4 iFAIL
C      
      INTEGER*4 ii, jj, kk
      INTEGER*4 NVBOOL
      PARAMETER (NVBOOL=100000000)
      INTEGER*4 vbool(1:NVBOOL)
C
      iFAIL = 0
C-----Check that supplied vector length for kprimes is allowable
      IF(kmax .GT. NVBOOL) THEN
          iFAIL = 1
          GOTO 80
      ENDIF
C-----Check that incoming kfac contains just integers of 2 or above...
      DO 10 jj = 1, nkfac
          IF(kfac(jj) .LT. 2) THEN
              iFAIL = 2
              GOTO 80
          ENDIF
10    CONTINUE
      DO 20 jj = 1, kmax
          vbool(jj) = 1
          kprimes(jj) = 0
20    CONTINUE
      kkmax = 0
      DO 32 jj = 1, nkfac
          ! integer divide drops the fractional part...
          DO 30 ii = 2, kmax/kfac(jj)
              kk = ii*kfac(jj)
              vbool(kk) = 0    ! kkth is not a prime
              IF( kk .GT. kkmax ) kkmax = kk
30        CONTINUE
32    CONTINUE
C
      jj = 0
      DO 40 ii = 2, kmax
          IF(vbool(ii) .EQ. 1) THEN
              jj = jj+1
              kprimes(jj) = ii     ! vector of primes
          ENDIF
40    CONTINUE
      nprimes = jj                 ! the number of primes up to kmax
C     
80    jj = jj
      END ! FINDP