Welcome, guest | Sign In | My Account | Store | Cart

Basic idea was to see the difference between different prime algorithms in time. Also they are not perfect the output shows that really higher numbers let grow the difference why I have separated this into functions to make it visible. I add this here because I have been missing this very often when I have been searching for algorithms.

  • The 'is_prime' is a well known way of checkin for a number being prime or not.
  • The sieve of Erastothenes is simply to strike out multiples of a given value; the primes will remain.
  • the function 'profile' is a decorator for functions measuring the execution time
  • Some information are in the comments of the code
Python, 250 lines
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
"""
 Topic: prime algorithms.
 Here is an example output of this script:

 Totally free. No warranty.

Python 3.1.2 (r312:79149, Mar 21 2010, 00:41:52) [MSC v.1500 32 bit (Intel)]

 [calc_slower]
 seeking primes between (0..100000) ...
 ...profiling of 'create_primes_1': tooks 6.974707 seconds
 ...profiling of 'create_primes_2': tooks 0.195253 seconds
 ...profiling of 'create_primes_3': tooks 0.090312 seconds
 ...profiling of 'create_primes_4': tooks 0.062241 seconds
 ...profiling of 'create_primes_5': tooks 0.046717 seconds
 ...profiling of 'create_primes_6': tooks 0.024859 seconds
 done - 9592 primes found.

 [calc_medium]
 seeking primes between (0..1000000) ...
 ...profiling of 'create_primes_2': tooks 4.662721 seconds
 ...profiling of 'create_primes_3': tooks 1.133545 seconds
 done - 78498 primes found.

 [calc_faster]
 seeking primes between (0..10000000) ...
 ...profiling of 'create_primes_3': tooks 13.872371 seconds
 ...profiling of 'create_primes_4': tooks 7.766646 seconds
 ...profiling of 'create_primes_5': tooks 6.069688 seconds
 done - 664579 primes found.

Are there no better solutions? ... smile ... Of course there are many!
But unfortunately some of them require a much better understanding in
math and I'm afraid (for myself) that some explanations are not simple
enough so I can present them here.
"""
import sys
import time
import math

def profile(function):
    """ intended to be used as decorator for a function or a method to check
        the execution time """

    def decorate(args):
        """ concrete decorator measuring the execution time of
            given function or method """
        start = time.clock()
        result = function(args)
        print("...profiling of '%s': tooks %f seconds" \
              % (function.__name__, time.clock()-start))
        #print("  -> %s" % (function.__doc__))
        return result
    return decorate


def is_prime(val):
    """ simple check for a prime value.
        Note: as a function called several time it is expensive. """
    if val <= 1:
        return False
    if val == 2:
        return True
    if val % 2 == 0:
        return False

    for divisor in range(3, int(math.sqrt(val))+1, 2):
        if val % divisor == 0:
            return False

    return True


@profile
def create_primes_1(max_n):
    """ creating primes up to a maximum value. The even numbers are not
        touched. The primes itself are used to check for possible division. """
    primes = []
    for val in range(3, max_n+1, 2):
        found = False
        for divisor in primes:
            if val % divisor == 0:
                found = True
                break

        if not found:
            primes.append(val)
    return [2] + primes


@profile
def create_primes_2(max_n):
    """ creating primes up to a max. value calling 'is_prime' for each value. """
    primes = [2]+[val for val in range(3, max_n+1, 2) if is_prime(val)]
    return primes

@profile
def create_primes_3(max_n):
    """ using sieve of sundaram:
        http://en.wikipedia.org/wiki/Sieve_of_Sundaram
        The description was really understandable   """
    limit    = max_n + 1
    sieve    = [True] * (limit)
    sieve[1] = False

    for i in range(1, limit//2):
        f = 2 * i + 1
        for j in range(i, limit//2):
            k = i + j * f
            if k <= max_n:
                sieve[k] = False
            else:
                break

    return [2,3]+[2*k+1 for k in range(1, int(max_n/2)) if sieve[k]]

@profile
def create_primes_4(max_n):
    """ creating primes up to a maximum value using a sieve algorithm. All
        multiples of a prime are flagged as 'no prime'. In addition there
        is an optimization by ignoring values flagged as none prime when
        proceeding to next value."""
    sieve = [True] * (max_n+1) # default: all are primes
    sieve[0] = False           # init:    0 is not a prime
    sieve[1] = False           # init:    1 is not a prime

    val = 2
    while val <= max_n:
        factor = val
        # strike out values not being a prime
        noprime = val * factor
        while noprime <= max_n:
            sieve[noprime] = False
            factor += 1
            noprime = val * factor
        # next value
        val += 1
        while val <= max_n and sieve[val] == False:
            val += 1

    return [n for n in range(max_n) if sieve[n] == True]


@profile
def create_primes_5(max_n):
    """ creating primes up to a maximum value using a sieve algorithm. All
        multiples of a prime are flagged as 'no prime'. In addition there
        is an optimization by ignoring values flagged as none prime when
        proceeding to next value."""
    sieve  = [False, True] * (int(max_n / 2))
    sieve += [True]

    if not max_n % 2 == 0:
        sieve += [False]


    sieve[1] = False
    sieve[2] = True
    primes   = [2]

    val = 3
    while val <= max_n:
        # now we have one prime
        primes.append(val)
        # strike out values not being a prime
        noprime = val + val
        while noprime <= max_n:
            sieve[noprime] = False
            noprime += val
        # next value
        val += 1
        while val <= max_n and sieve[val] == False:
            val += 1

    return primes


@profile
def create_primes_6(max_n):
    """ creating primes up to a maximum value using a sieve algorithm. All
        multiples of a prime are flagged as 'no prime'. In addition there
        is an optimization by ignoring values flagged as none prime when
        proceeding to next value."""
    sieve  = [False, True] * (max_n // 2)
    sieve += [False]

    sieve[1] = False
    sieve[2] = True
    primes   = [2]

    val = 3
    while val <= max_n:
        # now we have one prime
        primes.append(val)
        # strike out values not being a prime
        offset  = val * 2
        noprime = val + offset
        while noprime <= max_n:
            sieve[noprime] = False
            noprime += offset
        # next value
        val += 2
        while val <= max_n:
            if not sieve[val]:
                val += 2
            else:
                break

    return primes


def calc_slower(max_n):
    """ scenario I - displaying relation and validate results """
    print("\n[%s]" % (sys._getframe().f_code.co_name))
    print("seeking primes between (0..%d) ..." % (max_n))
    primes1 = create_primes_1(max_n)
    primes2 = create_primes_2(max_n)
    primes3 = create_primes_3(max_n)
    primes4 = create_primes_4(max_n)
    primes5 = create_primes_5(max_n)
    primes6 = create_primes_6(max_n)

    assert primes1 == primes2 == primes3 == primes4 == primes5 == primes6
    print("done - %d primes found." % (len(primes4)))


def calc_medium(max_n):
    """ scenario II - testing 'is_prime' alone """
    print("\n[%s]" % (sys._getframe().f_code.co_name))
    print("seeking primes between (0..%d) ..." % (max_n))
    primes2 = create_primes_2(max_n)
    primes3 = create_primes_3(max_n)
    print("done - %d primes found." % (len(primes2)))


def calc_faster(max_n):
    """ scenario III - testing sieve alone """
    print("\n[%s]" % (sys._getframe().f_code.co_name))
    print("seeking primes between (0..%d) ..." % (max_n))
    primes3 = create_primes_3(max_n)
    primes4 = create_primes_4(max_n)
    primes5 = create_primes_5(max_n)
    print("done - %d primes found." % (len(primes4)))

if __name__ == "__main__":
    print("Python %s" % (sys.version))
    max_n = 100000
    calc_slower(max_n      )
    calc_medium(max_n *  10)
    calc_faster(max_n * 100)

Python is a really nice language and I like to see that you can calculate expensive things with an interpreter. One goal why I have been providing this here is that I would like to hear the improvements but - please - do not just provide code. Give an explanation; I'm willing to learn. I have seen other sieves but it does not make sense just to copy code if you (I) can not understand.

I know from the sieve that there is a variant with squaring but - please - can someone explain this a little bit?

Interesting links I have found:

If you can add something, please tell.

2 comments

Thomas Lehmann (author) 13 years, 9 months ago  # | flag

I have changed the create_primes_5 function to be faster. One issue was to avoid touching the multiple of 2's (see: offset). Another issue was the last loop which tries to find the next prime - even numbers can be ignored so I have changed the step from 1 to 2.

Thomas Lehmann (author) 13 years, 9 months ago  # | flag

New sieve added: sieve of Sundaram I've written the code myself using this page: http://en.wikipedia.org/wiki/Sieve_of_Sundaram The explanation there has been quite well. The optimization I did are not perfect; I'm not quite sure about the outer limits to avoid to much calculation. At least the assertions ensure that the result is the same as for other calculations.