Department of Mathematics, University of Utah
	Department of Mathematics, College of Science, University of Utah, Salt Lake City, Utah 84112-0090, USA. 801-581-6851, 801-581-4148 (fax)

A FORTRAN code

You can cut and paste the FORTRAN code below to count prime numbers. Otherwise this page is quite boring.

On most machines the program will run in 32 bit arithmetic, letting you analyze numbers up to

You may be able to run the program with 64 bit arithmetic on your machine after replacing the implicit integer statements in the program with corresponding implicit integer*8 statements. In that case you can run the program for values of n up to

      implicit integer (a - z)
      common             primes,      twins,       goldbach
      parameter (maxsize = 1000000)
      logical            top(maxsize +1),          bottom(maxsize +1)
      primes = 1
      twins = 0
      goldbach = 0
      write (*, *)' give an even integer greater than 4'                -
      read (*, *) n
      if (mod(n, 2) .eq. 1) n = n + 1
      if (n .lt. 4) n = 4
      twin1 = 0
      twin2 = 0
      if (n .le. 2*maxsize) then
         n1 = 3
         call sieve(n1, n, bottom, p1, p2)
         do 10 k = 1, n/4
            if (bottom(k) .and. bottom(n/2 - k - 1)) then
               goldbach = goldbach + 1
               g1 = 2*k + 1
               g2 = n - 2*k - 1
            end if
   10    continue
         lngth = (n - n1)/2
         do 20 i = lngth, 2, - 1
            if (bottom(i) .and. bottom(i - 1)) then
               twin1 = n1 + 2*(i - 2)
               twin2 = twin1 + 2
               go to 30
            end if
   20    continue
   30    continue
      else
         steps = (n - 2)/(4*maxsize)
         do 70 i = 1, steps
            n1 = 2*(i - 1)*maxsize
            n2 = n1 + 2*maxsize
            nlast = n2
            call sieve(n1, n2, bottom, p1, p2)
            if (i .gt. 1) then
               if (p1 - lasttop .eq. 2) then
                  twins = twins + 1
               end if
            end if
            lasttop = p2
            nn = n1
            n1 = n - n2
            n2 = n - nn
            if (i .eq. 1) n2 = n
            call sieve(n1, n2, top, p1, p2)
            if (i .eq. 1) then
               do 40 k = maxsize, 2, - 1
                  if (top(k) .and. top(k - 1)) then
                     twin1 = n1 + 2*(k - 2) + 1
                     twin2 = twin1 + 2
                     go to 50
                  end if
   40          continue
   50          continue
            end if
            if (i .gt. 1) then
               if (lastbot - p2 .eq. 2) then
                  twins = twins + 1
               end if
            end if
            lastbot = p1
            do 60 k = 1, maxsize
               if (bottom(k) .and. top(maxsize - k + 1)) then
                  goldbach = goldbach + 1
                  g1 = nn + 2*k - 1
                  g2 = n2 - 2*k + 1
               end if
   60       continue
   70    continue
         call sieve(nlast, n1, top, p1, p2)
      end if
      do 80 k = 1, (n1 - nlast)/4
         if (top(k) .and. top((n1 - nlast)/2 - k + 1)) then
            goldbach = goldbach + 1
            g1 = nlast + 2*k - 1
            g2 = n1 - 2*k + 1
         end if
   80 continue
      write (*, 1000) n, primes
      write (*, 1010) n, twins
      write (*, 1020) n, twin1, twin2
      write (*, 1030) n, goldbach, n, g1, g2
      stop
 1000 format (/' The number of primes <= ', i12, ' is ', i12)
 1010 format (/' The number of twins below ', i12, ' is ', i12)
 1020 format (' The prime twin closest to but below ', i12, ' is ',
     +   2i12)
 1030 format (/i12, ' can be written in ', i12, ' different ways as', /
     +   ' the sum of two odd primes, for example,'/, i12, '=', i12,
     +   '+', i12)
      end
*
*
      subroutine  sieve(n1, n2, l, p1, p2)
      implicit integer (a - z)
      common             primes,      twins,       goldbach
      integer            n1,          n2
      logical            l(1),        first
*
* make bounds odd
*
      nn1 = n1 + 1 - mod(n1, 2)
      nn2 = n2 - 1 + mod(n2, 2)
*
* initialize l
*
      lngth = (nn2 - nn1)/2 + 1
      do 10 i = 1, lngth
         l(i) = .true.
   10 continue
      if (n1 .eq. 0 .or. n1 .eq. 1) l(1) = .false.
*
* get ready for deleting non-primes
*
      nsqrt = nint(sqrt(float(n2)))
      do 30 i = 3, nsqrt, 2
         if (mod(i, 3) .eq. 0 .and. i .gt. 3) go to 30
         if (mod(i, 5) .eq. 0 .and. i .gt. 5) go to 30
         if (mod(i, 7) .eq. 0 .and. i .gt. 7) go to 30
         if (mod(i, 11) .eq. 0 .and. i .gt. 11) go to 30
         if (mod(i, 13) .eq. 0 .and. i .gt. 13) go to 30
*
* calculate beginning index
*
         ifirst = nn1 - mod(nn1, i)
         if (ifirst .lt. nn1) ifirst = ifirst + i
         if (mod(ifirst, 2) .eq. 0) ifirst = ifirst + i
         ifirst = (ifirst - nn1)/2 + 1
         if (i .eq. nn1 + (ifirst - 1)*2) ifirst = ifirst + i
         do 20 j = ifirst, lngth, i
            l(j) = .false.
   20    continue
   30 continue
*
* actually count primes
*
      first = .false.
      do 40 j = 1, lngth
         if (l(j)) then
            primes = primes + 1
            if ( .not. first) then
               p1 = nn1 + 2*(j - 1)
               first = .true.
            end if
            p2 = nn1 + 2*(j - 1)
         end if
   40 continue
*
* count twins in this array
*
      do 50 i = 2, lngth
         if (l(i) .and. l(i - 1)) then
            twins = twins + 1
         end if
   50 continue
      return
      end