|
|
|
Eratosthenes was a Greek Mathematician, who lived 276 BC to 194 BC in northern Africa. He studied prime numbers, and is credited with the invention of the method of finding them illustrated on this page, called the Seive of Eratosthenes.
In this method, numbers are written in a table, such as the one, below. Beginning with the first prime, 2, all "proper" multiples of 2 -- that is, multiples other than 2 itself -- are crossed off. The next number not already crossed off is guaranteed to be prime -- 3 in this case. Continue crossing off the proper multiples of each number. Since this table is square, when you have finished crossing off all the proper multiples of the primes in the first row, the table will contain only primes.
Start with the number 2.
When all proper multiples of all numbers in the first row are deleted, the table will contain only primes
Do you understand what a prime number is? If not, send me an email (click the link, below) and I will be glad to help you. If so, then the next topic is finding the Prime Factorization of a number.
If you need to process a whole bunch of primes (e.g. if you want to count them, or count the ones with certain properties, etc.) then an efficient sieve is the way to go. One of the ways to improve the efficiency is to reduce the storage needed for the array itself. Of course, besides 2, only odd numbers can be primes, cutting out half the storage. In addition, besides 3, multiples of 3 can't be prime. These facts result in an efficient way to reduce the size of the sieve by a factor of 3.
The next step in efficiency comes from the fact that, apart from the first three primes 2, 3, and 5, in every group of 60 numbers, only the ones equivalent (mod 60) to 1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, or 59. It's handy that there are 16 possible primes in every group of 60, because then the sieve can use a single 16-bit integer to represent each such group.
The first part of the program runs the sieve, loading up the bits of the integers in the array p(), which is the sieve. Since p() starts out all zeros, we actually set the bits for the non-primes. After the sieve has been loaded, each element of p represents a group of 60 possible prime numbers. So to find out if i is prime, we split it into i1 (group number) and i2 (number within group). If fb(i2)≠0, then i2 is a possible prime (one of 1, 7, 11, 13, 17, etc.) and so then it's worth checking whether the bit of p(i1) for this prime is still zero by calculating the logical "and" of p(i1) with fb(i2).
Code to check whether i is prime:
i1 = Int(i / 60)
i2 = i - CDbl(i1) * 60
If fb(i2) <> 0 And (fb(i2) And p(i1)) = 0 Then countPrimes = countPrimes + 1
Entire program to check how many ways a number, nmax, can be represented as the sum of two primes:
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
www.faust.fr.bw.schule.de/mhb/eratosiv.htm is the source of the JavaScript used to make this page.
www-history.mcs.st-andrews.ac.uk/history/Mathematicians/Eratosthenes.html gives the history of Eratosthenes.
How to find the Prime Factorization of a number.
The webmaster and author of this Math Help site is Graeme McRae.