Source code in the fortran programming language

      DO WHILE(F*F <= LST)         !But, F*F might overflow the integer limit so instead,
      DO WHILE(F <= LST/F)                      !Except, LST might also overflow the integer limit, so 
      DO WHILE(F <= (IST + 2*(SBITS - 1))/F)    !Which becomes...
      DO WHILE(F <= IST/F + (MOD(IST,F) + 2*(SBITS - 1))/F) !Preserving the remainder from IST/F.


      MODULE PRIMEBAG	!Need prime numbers? Plenty are available.
C   Creates and expands a disc file for a sieve of Eratoshenes, representing odd numbers only and starting with three.
C   Storage requirements: an array of N prime numbers in 16/32/64 bits vs. a bit array up to the 16/32/64 bit limit.
C Word size               N            Prime   N words in bits           Bit array in bits.
C     8 bit            P(31) =           127               248                         128
C                      P(54) =           251               432                         256
C    16 bit         P(3,512) =        32,749            56,192                      32,768
C                   P(6,542) =        65,521           104,672                      65,536
C    32 bit   P(105,097,565) = 2,147,483,647     3,363,122,080               2,147,483,648
C             P(203,280,221) = 4,294,967,291     6,504,967,072               4,294,967,296
C    64 bit        2.112E17                ?          1.352E19   9,223,372,036,854,775,808 ~ 9.22E18
C from n/Ln(n)     4.158E17                ?          2.661E19  18,446,744,073,709,551,616 ~ 1.84E19
       INTEGER MSG	!I/O unit number.
       INTEGER SSTASH	!For attachment to my stash file.
       INTEGER SRECLEN,SCHARS,SBITS	!Sizes.
       INTEGER SORG	!Where the sieve starts. This must be three.
       INTEGER SLAST	!Last record in my stash file.
       DATA SSTASH,SREC,SLAST/0,0,0/	!Prepared by PRIMEBAG.
       PARAMETER (SRECLEN = 1024)	!4K disc bloc size, but RECL (in OPEN) is in terms of four-byte integers.
       PARAMETER (SCHARS = (SRECLEN - 1)*4)	!Reserving space for one number at the start.
       PARAMETER (SBITS = SCHARS*8)	!Known size of a character.
       PARAMETER (SORG = 3)		!First odd number past two, which is not odd.
       CHARACTER*(*) SFILE		!A name is needed.
       PARAMETER (SFILE = "C:/Nicky/RosettaCode/Primes/PrimeSieve.bit")	!I don't have to count the characters.
Components of a buffered record for the stash.
       INTEGER SREC	!The record number.
       CHARACTER*1 C4(4)	!The start of the record - a counter.
       CHARACTER*1 SCHAR(0:SCHARS - 1)	!The majority of the record - a bit array, packed in 8-bit blobs...
Collect some bit twiddling assistants for AND and OR, rather than bit shifting.
       CHARACTER*1 BITON(0:7),BITOFF(0:7)	!Functions IBSET and IBCLR may not be available, and are little-endian anyway.
       PARAMETER (BITON =(/CHAR(2#10000000),CHAR(2#01000000),	!128,  64,	Reading strictly left-to-right.
     1                     CHAR(2#00100000),CHAR(2#00010000),	! 32,  16,	Uncompromising bigendery.
     1                     CHAR(2#00001000),CHAR(2#00000100),	!  8,   4,	Not just for bytes in words,
     3                     CHAR(2#00000010),CHAR(2#00000001)/))	!  2,   1.	But also bits in bytes.
       PARAMETER (BITOFF=(/CHAR(2#01111111),CHAR(2#10111111),	!127, 191,	BITON + BITOFF = 255.
     2                     CHAR(2#11011111),CHAR(2#11101111),	!223, 239,
     1                     CHAR(2#11110111),CHAR(2#11111011),	!247, 251,
     3                     CHAR(2#11111101),CHAR(2#11111110)/))	!253, 254.
       CONTAINS
        INTEGER FUNCTION I4UNPACK(C4)	!Convert four successive characters into an integer.
         CHARACTER*1 C4(4)	!The characters.
          I4UNPACK = ((ICHAR(C4(1))*256 + ICHAR(C4(2)))*256	!Convert the first four bytes
     1               + ICHAR(C4(3)))*256 + ICHAR(C4(4))		!To a four-byte integer.
        END FUNCTION I4UNPACK	!Big-endian style, irrespective of cpu endianness.
        SUBROUTINE C4PACK(I4)	!Convert an integer into successive bytes.
Could return the result via a fancy function, but for now a global variable will do.
         INTEGER I4,N	!The integer, and a copy to damage.
         INTEGER I	!A stepper.
          N = I4	!Keep the original safe.
          DO I = 4,1,-1	!Know that four characters will do. Fixed format makes this easy.
            C4(I) = CHAR(MOD(N,256))	!Grab the low-order eight bits.
            N = N/256			!And shift right eight.
          END DO			!Do it again.
        END SUBROUTINE C4PACK	!Stored big-endianly, irrespective of cpu endianness.

        LOGICAL FUNCTION GRASPPRIMEBAG(F)
         INTEGER F		!The I/O unit number to use.
         LOGICAL EXIST		!Use the keyword as a name
         INTEGER IOSTAT		!And don't worry over assignment direction.
         CHARACTER*3 STYLE	!One way or another.
          SSTASH = F		!I shall use it.
          INQUIRE (FILE = SFILE,EXIST = EXIST)	!Trouble with a missing "path" may arise.
          IF (EXIST) THEN	!If the file exists,
            STYLE = "OLD"	!I shall read it.
           ELSE			!But if it doesn't,
            STYLE = "NEW"	!I shall create it.
          END IF		!Enough prevarication.
          OPEN(SSTASH,FILE = SFILE, STATUS = STYLE,	!Go for the file.
     &     ACCESS = "DIRECT", RECL = SRECLEN, FORM = "UNFORMATTED",	!I have plans.
     &     ERR = 666, IOSTAT = IOSTAT)					!Which may be thwarted.
          IF (EXIST) THEN	!If there is one...
            CALL READSCHAR(1)	!The first record is also a header.
            SLAST = I4UNPACK(C4)	!The number of records stored.
           ELSE			!Otherwise, start from scratch.
            SLAST = 0			!No saved records.
            CALL PSURGE(SCHAR)		!During preparation of the first batch of bits.
          END IF		!All should now be in readiness.
          GRASPPRIMEBAG = .TRUE.!So, feel confidence.
         RETURN			!And escape.
  666     WRITE (*,667) IOSTAT,SFILE	!But, something may have gone wrong.
  667     FORMAT ("Pox! Error code ",I0,	!A "hole" in the directory path?
     1     " when attempting to open file ",A)	!Read-only access allowed when I want "update"?
          GRASPPRIMEBAG = .FALSE.		!Whatever, it didn't work.
        END FUNCTION GRASPPRIMEBAG	!So much for that.

        SUBROUTINE READSCHAR(R)	!Get record R into SCHAR, which may already hold it.
         INTEGER R		!The record number desired.
          IF (R.EQ.SREC) RETURN	!Perhaps it is already to hand.
          SREC = R		!If not, move attention to it.
          READ (SSTASH,REC = SREC) C4,SCHAR	!And read the record.
        END SUBROUTINE READSCHAR!Thus, I have a buffer too.

        LOGICAL FUNCTION PSURGE(BIT8)	!Add another record to the stash.
C   Surges forward into the next batch of primes, to be stored via a bit array in the file.
C   Each record starts with a count of the number of primes that have gone before.
C   Except that for the first record, this is the record counter for the stash file.
C   Except that when starting the second record, one is also the number of primes before SORG.
         CHARACTER*1 BIT8(0:SCHARS - 1)	!Watch out! This may be SCHAR itself!
         INTEGER IST,LST	!The numbers spanned by the surge.
         INTEGER F		!A factor.
         INTEGER I		!Another factor and a stepper.
         INTEGER C		!Index for array BIT8.
         INTEGER NP		!Number of primes.
Carry forward the count of previous primes to start the following record..
   10     IF (SLAST.GT.0) THEN	!Is there a previous record?
            CALL READSCHAR(SLAST)	!Yes. Grab it. A good chance this is already in C4,SCHAR.
            NP = I4UNPACK(C4)		!Its count of the primes accumulated before it.
            DO I = 0,SCHARS - 1		!Find out how namy primes it fingered by scanning its bits.
              NP = NP + COUNT(IAND(ICHAR(SCHAR(I)),ICHAR(BITON)).NE.0)	!Whee! Eight at a go!
            END DO			!On to the next byte.
          END IF		!When creating a new record, its follower may not be sought in this run.
Concoct the next batch of bits. Contorted calculations avoid integer overflow.
   20     BIT8 = CHAR(255)		!All bits are aligned with numbers that might prove to be prime.
          IST = SORG + SLAST*(2*SBITS)	!Bit(0) of BIT8(0) corresponds to IST.
          LST = IST + 2*(SBITS - 1)	!Bit(last) to this number. Remember, only odd numbers have bits.
          IF (IST.LE.0) THEN		!Humm. I'd better check.
            WRITE (MSG,21) SLAST,IST,LST	!This works only with two's complement integers.
   21       FORMAT (/,"Integer overflow in the sieve of Eratosthenes!",	!Oh dear.
     1      /,"Advancing from surge ",I0," to span ",I0," to ",I0)	!These numbers will look odd.
            PSURGE = .FALSE.			!But it is better than no indication of what went wrong.
           RETURN			!Give in.
          END IF		!Enough worrying.
          F = 3			!The first possible factor. Zapping will start at F²
c         DO WHILE(F.LE.LST/F)	!If F² is past the end, so will be still larger F: enough.
          DO WHILE(F.LE.IST/F + (MOD(IST,F) + 2*(SBITS - 1))/F)	!"Synthetic division" avoiding overflow.
            I = (IST - 1)/F + 1		!I want the first multiple of F in IST:LST. F may be a factor of IST.
            IF (MOD(I,2).EQ.0) I = I + 1!If even, advance to the next odd multiple. Even numbers are omitted by design.
            IF (I.LT.F) I = F		!Less than F is superfluous: the position was zapped by earlier action.
c           I = (I*F - IST)/2			!Current bit positions are for IST, IST+2, IST+4, etc.
            I = ((I - IST/F)*F - MOD(IST,F))/2	!Avoids overflow when calculating the start value, I*F.
            DO I = I,SBITS - 1,F	!Zap every F'th bit along. This is the sieve of Eratosthenes.
              C = I/8				!Eight bits per character.
              BIT8(C) = CHAR(IAND(ICHAR(BIT8(C)),	!For F = 3 and 5, characters will be hit more than once.
     1                            ICHAR(BITOFF(MOD(I,8)))))	!Whack a bit. All the above just for this!
            END DO			!On to the next bit.
   22       F = NEXTPRIME(F)		!So much for F. Next, please.
          END DO		!Are we there yet?
Correct the count in the header, if this is an added record.
   30     IF (SLAST.GT.0) THEN	!So, was there a pre-existing header record?
            CALL READSCHAR(1)		!Yes. Get the header record into C4,SCHAR.
            CALL C4PACK(SLAST + 1)	!This is the new record count.
            WRITE (SSTASH,REC = 1) C4,SCHAR	!Write it all back.
            SCHAR = BIT8	!Ensure that SCHAR and SREC will be agreed.
          END IF		!So much for the header's count.
Cast the bits into the stash by writing record SLAST + 1..
   40     IF (SLAST.EQ.0) THEN	!If we're writing the first record,
            CALL C4PACK(1)		!Then this is the record count.
           ELSE			!Otherwise,
            CALL C4PACK(NP)		!Place the previous primes count.
          END IF		!All this to help PRIME(i).
          SLAST = SLAST + 1	!This is now the last stashed record.
          WRITE (SSTASH,REC = SLAST) C4,BIT8	!I/O directly from the work area?
          SREC = SLAST		!This is where BIT8 was written.
          PSURGE = .TRUE.	!That assumes BIT8 is not SCHAR for SLAST > 1.
        END FUNCTION PSURGE	!That was fun!

        RECURSIVE SUBROUTINE GETSREC(R)	!Make present the bit array belonging to record R.
         INTEGER R		!The record number..
         CHARACTER*1 BIT8(0:SCHARS - 1)	!A scratchpad. Others may be relying on SCHAR.
          IF (SLAST.LE.0) RETURN!DANGER! The first record is being initialised!
          DO WHILE(SLAST.LT.R)	!If we haven't reached so far,
            IF (.NOT.PSURGE(BIT8)) THEN	!Slog forwards one record's worth.
              WRITE (MSG,1) R			!Or maybe not.
    1         FORMAT ("Cannot prepare surge ",I0)	!Explain.
              STOP "No bits, no go."			!And quit.
            END IF			!And having prepared the next block of bits,
          END DO		!Check afresh.
          CALL READSCHAR(R)	!Read the desired record's bits.
        END SUBROUTINE GETSREC	!Done.

        INTEGER FUNCTION PRIME(N)	!P(1) = 2, P(2) = 3, etc.
C   Calculate P(n) ~ n.ln(n)
C                  ~ n{ln(n) + ln(ln(n)) - 1 + (ln(ln(n)) - 2)/ln(n) - [ln(ln(n))**2 - 6*log(log(n)) + 11]/[2*(ln(n))**2] + ....}
C   J.B.Rosser's 1938 Theorem: n[ln(n) + ln(ln(n)) - 1] < P(n) < n[ln(n) + ln(ln(n))]
C    or, with E = ln(n) + ln(ln(n)),           n[E - 1] < P(n) < n[E]
C   Experimentation shows that the undershoot of the first two terms involves many records worth of bits.
C   Including additional terms does much better, but can overshoot.
         INTEGER N		!The desired one.
         INTEGER R,NP		!Counts.
         INTEGER B,C		!Bit and character indices.
         DOUBLE PRECISION EST,LN,LLN	!Hope, if not actuality.
          IF (N.LE.0) STOP "Primes are counted positively!"	!Something must be wrong!
          IF (N.LE.1) THEN	!The start of the bit array being preempted.
            PRIME = 2			!So, no array access.
           ELSE			!Otherwise, the fun begins.
            LN = LOG(DFLOAT(N))	!Here we go.
            LLN = LOG(LN)	!A popular term.
            EST = N*(LN		!Estimate the value of the N'th prime.
     1               + LLN - 1		!Second term
     2               + (LLN - 2)/LN		!Third term.
     3               - (LLN**2 - 6*LLN + 11)/(2*LN**2))	!Fourth term.
            R = (EST - SORG)/(2*SBITS) + 1	!Thereby selecting a record to scan.
            IF (R.LE.0) R = 1			!And not making a mess with N < 6 or so.
    9       CALL GETSREC(R)	!Go for the record.
            IF (R.LE.1) THEN	!The first record starts with the record count.
              NP = 1		!And I know how many primes precede its start point
             ELSE		!While for all subsequent records,
              NP = I4UNPACK(C4)	!This counts the number of primes that precede record R's start number.
            END IF		!So now I'm ready to count onwards.
            IF (N.LE.NP) THEN	!Maybe not.
              R = R - 1			!The estimate took me too far ahead.
              GO TO 9			!Try again.
            END IF		!Could escalate to a binary search or even an interpolating search.
Commence scanning the bits.
            C = 0		!Start with the first character of SREC..
            B = -1		!Syncopation. The formula is known to always under-estimate.
   10       IF (NP.LT.N) THEN	!Are we there yet?
   11         B = B + 1			!No. Advance to the next bit.
              IF (B.GE.8) THEN		!Overflowed a character yet?
                B = 0			!Yes. Start afresh at the first bit.
                C = C + 1		!And advance one character.
                IF (C.GE.SCHARS) THEN	!Overflowed the record yet?
                  C = 0			!Yes. Start afresh at its first character.
                  R = R + 1		!And advance to the next record.
                  CALL GETSREC(R)	!Possibly, create it.
                END IF			!So much for records.
              END IF			!We're now ready to test bit B of character C of record R.
              IF (IAND(ICHAR(SCHAR(C)),ICHAR(BITON(B))).EQ.0) GO TO 11	!Not a prime. Search on.
              NP = NP + 1		!Count another prime.
              GO TO 10			!Pehaps this will be the one.
            END IF		!So much for the search.
            PRIME = SORG + (R - 1)*(2*SBITS) + (C*8 + B)*2	!The corresponding number.
            IF (PRIME.LE.0) WRITE (MSG,666) N,PRIME	!Or, possibly not.
  666       FORMAT ("Integer overflow! Prime(",I0,") gives ",I0,"!")	!Let us hope the caller notices.
          END IF		!So, all going well,
        END FUNCTION PRIME	!It is found.

        RECURSIVE INTEGER FUNCTION NEXTPRIME(N)	!Keep right on to the end of the road.
Can invoke GETSREC, which can invoke PSURGE, which ... invokes NEXTPRIME. Oh dear.
         INTEGER N	!Not necessarily itself a prime number.
         INTEGER NN	!A value to work with.
         INTEGER R	!A record number into the stash.
         INTEGER I,IST	!Number offsets.
         INTEGER C,B	!Character and bit index.
          IF (N.LE.1) THEN	!Suspicion prevails.
            NN = 2			!This is not represented in my bit array.
           ELSE			!Otherwise, the fun begins.
            NN = N + 1			!Advance, with a copy I can mess with.
            IF (MOD(NN,2).EQ.0) NN = NN + 1	!Thus, NN is now odd.
            IF (NN.LE.0) GO TO 666		!But perhaps not proper, due to overflow.
            R = (NN - SORG)/(2*SBITS)		!SORG is odd, so (NN - SORG) is even.
            CALL GETSREC(R + 1)		!The first record is numbered one, not zero.
            IST = SORG + R*(2*SBITS)	!The number for its first bit: even numbers are omitted..
            I = (NN - IST)/2		!Offset into the record. NN - IST is even.
            C = I/8			!Which character in SCHAR(0:SCHARS - 1)?
            B = MOD(I,8)		!Which bit in SCHAR(C)?
   10       IF (IAND(ICHAR(SCHAR(C)),ICHAR(BITON(B))).EQ.0) THEN	!On for a prime.
              NN = NN + 2	!Alas, it is off, so NN is not a prime. Perhaps this will be.
              B = B + 1		!Advance one bit. Each bit steps two.
              IF (B.GE.8) THEN	!Past the end of the character?
                B = 0			!Yes. Back to bit zero.
                C = C + 1		!And advance one chracter.
                IF (C.GE.SCHARS) THEN	!Past the end of the record?
                  IF (NN.LE.0) GO TO 666!Yes. If NN has overflowed, the end of the rope is reached.
                  C = 0			!Back to the start of a record.
                  R = R + 1		!Advance one record.
                  CALL GETSREC(R + 1)	!And read it. (Count is from 1, not 0).
                END IF		!So much for overflowing a record.
              END IF		!So much for overflowing a character.
              GO TO 10	!Try again.
            END IF		!So much for the bit array.
          END IF		!If there had been a scan.
          NEXTPRIME = NN	!The number for which the scan stopped.
          IF (NN.GT.0) RETURN	!All is well.
  666     WRITE (MSG,667) N,NN	!Or, maybe not. Careful: this won't appear if NEXTPRIME is invoked in a WRITE list.
  667     FORMAT ("Integer overflow! NextPrime(",I0,") gives ",I0,"!")	!The recipient could do a two's complement.
          NEXTPRIME = NN	!Prefer to return the bad value rather than fail to return anything.
        END FUNCTION NEXTPRIME	!No divisions, no sieving. Here, anyway

        INTEGER FUNCTION PREVIOUSPRIME(N)	!If N is good, this can't overflow.
         INTEGER N	!The number, not necessarily a prime.
         INTEGER NN	!A value to mess with.
         INTEGER R	!A record number.
         INTEGER I	!Offset.
         INTEGER C,B	!Character and bit fingers.
          IF (N.LE.3) THEN	!Suppress annoyances.
            NN = 2		!This is now called the first prime, not one.
           ELSE			!Otherwise, some work is to be done.
            NN = N - 1			!Step back one to ensure previousness.
            IF (MOD(NN,2).EQ.0) NN = NN - 1	!And here, oddness is a minimal requirement.
            R = (NN - SORG)/(2*SBITS)	!Finger the record containing the bit for NN.
            CALL GETSREC(R + 1)		!Record counting starts with one.
            I = (NN - (SORG + R*(2*SBITS)))/2	!Offset into that record.
            C = I/8			!Finger the character in SCHAR.
            B = MOD(I,8)		!And the bit within the character.
   10       IF (IAND(ICHAR(SCHAR(C)),ICHAR(BITON(B))).EQ.0) THEN	!On for a prime.
              NN = NN - 2	!Alas, it is off, so NN is not a prime. Perhaps this will be.
              B = B - 1		!Retreat one bit. Each bit steps two.
              IF (B.LT.0) THEN	!Past the start of the character?
                B = 7			!Yes. Back to the last bit.
                C = C - 1		!And retreat one chracter.
                IF (C.LT.0) THEN	!Past the start of the record?
                  C = SCHARS - 1	!Yes. Back to the end of a record.
                  R = R - 1		!Retreat one record.
                  CALL GETSREC(R + 1)	!And read it. (Count is from 1, not 0).
                END IF		!So much for overflowing a record.
              END IF		!So much for overflowing a character.
              GO TO 10		!Try again.
            END IF		!So much for the bit array.
          END IF		!Possibly, it was not needed.
          PREVIOUSPRIME = NN	!There.
        END FUNCTION PREVIOUSPRIME	!Doesn't overflow, either.

        LOGICAL FUNCTION ISPRIME(N)	!Could fool around explicity testing 2 and 3 and say 5,
         INTEGER N			!But that means also checking that N > 2, N > 3, and N > 5.
c        ISPRIME = N .EQ. NEXTPRIME(N - 1)	!This is so much easier, but involves scanning to reach the next prime.
         INTEGER R,IST,I,C,B		!Assistants for indexing the bit array.
          IF (N.LE.1) THEN	!First, preclude sillyness.
            ISPRIME = .FALSE.		!Not a prime.
          ELSE IF (N.EQ.2) THEN	!This is the only even number
            ISPRIME = .TRUE.		!That is a prime.
          ELSE IF (MOD(N,2).EQ.0) THEN	!Other even numbers
            ISPRIME = .FALSE.		!Are not prime numbers.
          ELSE			!Righto, now N is an odd number and there is a bit array for them.
            R = (N - SORG)/(2*SBITS)	!SORG is odd, so (N - SORG) is even.
            CALL GETSREC(R + 1)		!The first record is numbered one, not zero.
            IST = SORG + R*(2*SBITS)	!The number for its first bit: even numbers are omitted.
            I = (N - IST)/2		!Offset into the record. N - IST is even.
            C = I/8			!Which character in SCHAR(0:SCHARS - 1)?
            B = MOD(I,8)		!Which bit in SCHAR(C), indexing from zero?
            ISPRIME = IAND(ICHAR(SCHAR(C)),ICHAR(BITON(B))).GT.0	!The bit is on for a prime.
          END IF			!All that fuss to find a single bit.
        END FUNCTION ISPRIME		!But, no divisions up to SQRT(N) or the like.
      END MODULE PRIMEBAG	!Functions updating a disc file as a side effect...

      PROGRAM POKE
      USE PRIMEBAG
      INTEGER I,P,N,N1,N2	!Assorted assistants.
      INTEGER ORDER		!A collection of special values.
      PARAMETER (ORDER = 6)	!For one, two, and four byte integers.
      INTEGER EDGE(ORDER)	!Considered as two's complement and unsigned.
      PARAMETER (EDGE = (/31,54,3512,6542,105097565,203280221/))	!These primes are of interest.
      MSG = 6		!Standard output.

      IF (.NOT.GRASPPRIMEBAG(66)) STOP "Gan't grab my file!"	!Attempt in hope.

Case 1.
C     FORALL(I = 1:20) LIST(I) = PRIME(I) is rejected because function Prime(i) is rather impure.
   10 WRITE (MSG,11)
   11 FORMAT (19X,"First twenty primes: ", $)
      DO I = 1,20
        P = PRIME(I)
        WRITE (MSG,12) P
   12   FORMAT (I0,",",$)
      END DO

Case 2.
   20 WRITE (MSG,21)
   21 FORMAT (/,12X,"Primes between 100 and 150: ",$)
      P = 100
   22 P = NEXTPRIME(P)		!While (P:=NextPrime(P)) <= 150 do Print P;
      IF (P.LE.150) THEN	!But alas, no assignment within an expression.
        WRITE (MSG,23) P
   23   FORMAT (I0,",",$)
        GO TO 22
      END IF

Case 3.
   30 N1 = 7700	!Might as well parameterise this.
      N2 = 8000	!Rather than litter the source with explicit integers.
      N = 0
      P = N1
   31 P = NEXTPRIME(P)
      IF (P.LE.N2) THEN
        N = N + 1
        GO TO 31
      END IF
      WRITE (MSG,32) N1,N2,N
   32 FORMAT (/"Number of primes between ",I0," and ",I0,": ",I0)

Case 4.
   40 WRITE (MSG,41)
   41 FORMAT (/,"Tenfold steps...")
      N = 1
      DO I = 1,9	!This goes about as far as it can go.
        P = PRIME(N)
        WRITE (MSG,42) N,P
   42   FORMAT ("Prime(",I0,") = ",I0)
        N = N*10
      END DO

Cast forth some interesting values.
  100 WRITE (MSG,101)
  101 FORMAT (/,"Primes close to number sizes")
      DO N = 1,ORDER	!Step through the list.
        N1 = EDGE(N) - 1	!Syncopation for the special value.
        DO I = 1,2		!I want the prime on either side.
          N1 = N1 + 1			!So, there are two successive primes to finger.
          WRITE (MSG,102) N1			!Identify the index.
  102     FORMAT ("Prime(",I0,") = ",$)		!Piecemeal writing to the output,
          P = PRIME(N1)				!As this may fling forth a complaint.
          WRITE (MSG,103) P			!Show the value returned.
  103     FORMAT (I0,", ",$)			!Which may be unexpected.
        END DO			!On to the second.
        WRITE (MSG,*)		!End the line after the second result.
      END DO		!On to the next in the list.

      END	!Whee!


      P = NEXTPRIME(100)
      DO WHILE (P.LE.150)
        ...stuff...
        P = NEXTPRIME(P)
      END DO


 P:=100; WHILE (P:=NextPrime(P)) <= 150 DO stuff;

* incremental Sieve of Eratosthenes based on the paper,
* "Two Compact Incremental Prime Sieves"

      SUBROUTINE nextprime(no init, p)
      IMPLICIT NONE
      INTEGER*2, SAVE, ALLOCATABLE :: sieve(:,:)
      INTEGER, SAVE :: r, s, pos, n, f1, f2, sz
      INTEGER i, j, d, next, p, f3
      LOGICAL no init, is prime

      IF (no init) GO TO 10 
      IF (ALLOCATED(sieve)) DEALLOCATE(sieve)

* Each row in the sieve is a stack of 8 short integers. The
* stacks will never overflow since the product 2*3*5 ... *29
* (10 primes) exceeds a 32 bit integer. 2 is not stored in the sieve.

      ALLOCATE(sieve(8,3))
      sieve = reshape([(0_2, i = 1, 24)], shape(sieve))
      r = 3
      s = 9
      pos = 1
      sz = 1 ! sieve starts with size = 1
      f1 = 2 ! Fibonacci sequence for allocating new capacities
      f2 = 3 ! array starts with capacity 3
      n = 1
      p = 2  ! return our first prime
      RETURN

10    n = n + 2
      is prime = .true.
      IF (sieve(1, pos) .eq. 0) GO TO 20 ! n is non-smooth w.r.t sieve
      is prime = .false. ! element at sieve(pos) divides n
      DO 17, i = 1, 8

Clear the stack of divisors by moving them to the next multiple

         d = sieve(i, pos)
         IF (d .eq. 0) GO TO 20 ! stack is empty
         IF (d .lt. 0) d = d + 65536 ! correct storage overflow
         sieve(i, pos) = 0
         next = mod(pos + d - 1, sz) + 1

* Push divisor d on to the stack of the next multiple

         j = 1
12       IF (sieve(j, next) .eq. 0) GO TO 15
         j = j + 1
         GO TO 12
15       sieve(j, next) = d
17    CONTINUE

Check if n is square; if so, then add sieving prime and advance

20    IF (n .lt. s) GO TO 30
      IF (.not. is prime) GO TO 25
      is prime = .false. ! r = √s divides n
      next = mod(pos + r - 1, sz) + 1 ! however, r is prime, insert it.

      j = 1
22    IF (sieve(j, next) .eq. 0) GO TO 23
      j = j + 1
      GO TO 22
23    sieve(j, next) = r

25    r = r + 2
      s = r**2

Continue to the next array slot; grow the array by two when
* we get to the end to maintain the invariant size(sieve) > √n
* IF the size exceeds the array capacity, resize the arary.

30    pos = pos + 1
      IF (pos .le. sz) GO TO 40
      sz = sz + 2
      pos = 1

      IF (sz .le. f2) GO TO 40 ! so far, no need to grow
      f3 = f1 + f2
      f1 = f2
      f2 = f3
      sieve = reshape(sieve, [8, f2],
     &                pad = [(0_2, i = 1, 8*(f2 - f1))])

* Either return n back to the caller or circle back if n
* turned out to be composite.

40    IF (.not. is prime) GO TO 10
      p = n
      END SUBROUTINE


      INCLUDE 'sieve.f'

      PROGRAM RC Extensible Sieve
      IMPLICIT INTEGER (A-Z)

      WRITE (*, '(A)', advance='no')
     &   'The first 20 primes:'

      CALL nextprime(.false., p)
      DO 10, i = 1, 20
         WRITE (*, '(I3)', advance = 'no') p
10       CALL nextprime(.true., p)
      WRITE (*, *)

      WRITE (*, '(A)', advance = 'no')
     &   'The primes between 100 and 150:'

20    CALL nextprime(.true., p)
      IF (p .gt. 149) GO TO 30
      IF (p .gt. 99)
     &   WRITE (*, '(I4)', advance = 'no') p
      GO TO 20
30    WRITE (*, *)

      count = 0
40    CALL nextprime(.true., p)
      IF (p .gt. 7999) GO TO 50
      IF (p .gt. 7700) count = count + 1
      GO TO 40
50    WRITE (*, 100) count
100   FORMAT ('There are ', I0, ' primes between 7700 and 8000.')

      CALL nextprime(.false., p) ! re-initialize
      target = 1 ! target count
      n = 0 ! number of primes generated
60    n = n + 1
      IF (n .lt. target) GO TO 70
      WRITE (*, '(ES7.1,1X,I12)'), real(n), p
      IF (target .eq. 100 000 000) GO TO 80
      target = target * 10
70    CALL nextprime(.true., p)
      GO TO 60

80    END