Sub primeSieveOfEratosthenes()
   ' Bitmap "p" structure:
   '    If (p(n\60) And fb(n Mod 60)) = 0 then p is prime
   Dim nmax As Double ' size of the sieve; all primes up to nmax will be found.
   Dim nmaxsqrt As Long ' sqrt(nmax)
   Dim nmax60 As Long ' nmax\60
   Dim countPrimes As Long, countPairs As Long
   Dim p() As Integer ' sieve; each 16-bit integer has bits for a group of 60 possible primes
   Dim fb(59) As Integer ' "find bit" array, selects bit of p(n\60) that represents this prime; fb=0,1,0,0,0,0,0,2,...
   Dim ib(15) As Integer ' "index bit" array, indexes to the next nonzero element of fb; ib=1,7,11,13,...,53,59
   Dim fbib(59) As Integer ' fb(ib(i2)); fbib=1,2,4,8,...,-32768
   Dim i As Double, i1 As Long, i2 As Integer, j As Double, j1 As Long, j2 As Long
   Dim sec As Double ' number of seconds it takes me to update the sieve and count all the primes

   sec = Timer()
   nmax = Selection.Range("A1").Value ' use a billion for testing, or 179424673, or 2038074743
   nmaxsqrt = Int(Sqr(nmax))
   nmax60 = Int(nmax / 60)
   ReDim p(nmax60)
   ib(0) = 1: ib(1) = 7: ib(2) = 11: ib(3) = 13: ib(4) = 17: ib(5) = 19: ib(6) = 23: ib(7) = 29:
   ib(8) = 31: ib(9) = 37: ib(10) = 41: ib(11) = 43: ib(12) = 47: ib(13) = 49: ib(14) = 53: ib(15) = 59:
   fb(ib(0)) = 1: For i = 1 To 14: fb(ib(i)) = 2 * fb(ib(i - 1)): Next i: fb(ib(i)) = (-2) * fb(ib(i - 1))
   For i = 0 To 15: fbib(i) = fb(ib(i)): Next i

   i = 5: ' start off having counted 2, 3, and 5
   i1 = 0: i2 = 0
   Do
      ' Advance i = i1 * 60 + ib(i2) to the next prime
      Do
         For i2 = i2 + 1 To 15
            If (p(i1) And fbib(i2)) = 0 Then Exit Do
         Next i2
         i1 = i1 + 1
         If i1 > nmax60 Then GoTo finish
         i2 = -1
      Loop
      i = i1 * 60 + ib(i2)
      If i > nmaxsqrt Then GoTo finish
      ' knock out multiples of i, starting with i^2
      For j = i ^ 2 To nmax Step 2 * i
         j1 = Int(j / 60)
         'j2 = j - CDbl(j1) * 60
         p(j1) = p(j1) Or fb(j - CDbl(j1) * 60)
      Next j
   Loop
finish:
   countPrimes = 3 'start off having counted 2, 3, and 5
   For i = 5 To nmax Step 2
      i1 = Int(i / 60): i2 = i - CDbl(i1) * 60
      If fb(i2) <> 0 And (fb(i2) And p(i1)) = 0 Then countPrimes = countPrimes + 1
   Next i

   countPairs = 0
   i = 3
   j = nmax - i: j1 = Int(j / 60): j2 = j - CDbl(j1) * 60
   If fb(j2) <> 0 And (fb(j2) And p(j1)) = 0 Then countPairs = countPairs + 1
   For i = 5 To Int(nmax / 2) Step 2
      i1 = Int(i / 60): i2 = i - CDbl(i1) * 60
      If fb(i2) <> 0 And (fb(i2) And p(i1)) = 0 Then
         j = nmax - i: j1 = Int(j / 60): j2 = j - CDbl(j1) * 60
         If fb(j2) <> 0 And (fb(j2) And p(j1)) = 0 Then countPairs = countPairs + 1
      End If
   Next i

   Selection.Range("A1").Value = nmax
   Selection.Range("B1").Value = countPrimes
   Selection.Range("C1").Value = countPairs
   Selection.Range("D1").Value = Timer() - sec
   Selection.Range("A2").Select
End Sub