The sieve of Eratosthenes: you've been there; done that; have the T-shirt. The sieve of Eratosthenes was ancient history when Euclid was a schoolboy. You are ready for something less than 3000 years old. You are ready for The sieve of Sundaram.

Starting with the ordered set of +ve integers, mark every third starting at 4 (4;7;10...).

Step through the set and if the value is not marked output 2*n+1. So from 1 to 4 output 3 5 7.

4 is marked so skip for 5 and 6 output 11 and 13.

7 is marked, so no output but now also mark every fifth starting at 12 (12;17;22...)

as per to 10 and now mark every seventh starting at 17 (17;24;31....)

as per for every further third element (13;16;19...) mark every (9th;11th;13th;...) element.

The output will be the ordered set of odd primes.

Using your function find and output the first 100 and the millionth Sundaram prime.

The faithless amongst you may compare the results with those generated by The sieve of Eratosthenes.

References

 The article on Wikipedia.


 
      PROGRAM SUNDARAM
      IMPLICIT NONE
!
! Local variables
!
      INTEGER(8)  ::  curr_index
      INTEGER(8)  ::  i
      INTEGER(8)  ::  j
      INTEGER  ::  lim
      INTEGER(8)  ::  mid
      INTEGER  ::  primcount
      LOGICAL*1 , ALLOCATABLE , DIMENSION(:)  ::  primes !Array of booleans representing integers 
      lim = 10000000 ! Not the number of primes but the storage where the prime marker is held for the millionth prime
      ALLOCATE(primes(lim))
      primes(1:lim) = .TRUE.
                        !Set all to .True., we will later block out the known non-primes
      mid = lim/2
 
!Generate primes
      DO j = 1 , mid
         DO i = 1 , j
            curr_index = i + j + (2*i*j)
            IF( curr_index>lim )EXIT        ! Too big already, leave the loop.
            primes(curr_index) = .FALSE.    !This candidate will not produce a prime
         END DO
      END DO
!
      i = 0
      j = 0
      WRITE(6 , *)'The first 100 primes:'
      DO WHILE ( i < 100 )
         j = j + 1
         IF( primes(j) )THEN
            WRITE(6 , 34 , ADVANCE = 'no')j*2 + 1   !Take the candidate, multiply by 2, add 1, and you have a prime
 34         FORMAT(I0 , 1x)
            i = i + 1                       ! Counter used for printing
            IF( MOD(i,10)==0 )WRITE(6 , *)' '
         END IF
      END DO
! Now print the millionth prime
      primcount = 0
      DO i = 1 , lim
         IF( primes(i) )THEN
            primcount = primcount + 1
            IF( primcount==1000000 )THEN
               WRITE(6 , 35)'1 millionth Prime Found: ' , (i*2) + 1
 35            FORMAT(/ , a , i0)
               EXIT
            END IF
         END IF
      END DO
      DEALLOCATE(primes)
      STOP
      END PROGRAM SUNDARAM
 
 