The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.

The first thing that comes in mind is to use the previously defined `sieveEratosthenes`

function but where's the fun in that? So the very next thing that crossed my mind was to actually make that sluggish function a bit faster. If you do your research, you'll soon find a technique called wheel factorization.

Wheel factorization in short means that what you first do is that you pick a certain set of consecutive primes, then you multiply them together and in the range `1..result`

you identify multiples of your original primes (including the primes themselves). You're then supposed to write them around a circle but I don't see any good reason why not to write them in a line instead. Then in each new row (ring) the elements in marked spokes/columns can be discarded as multiples as well and only the remaining numbers are tested for primarility.

To give an example, this is what wheel factorization using 2 and 3 looks like:

` 1 `__2__ __3__ ~~4~~ 5 ~~6~~

` 7 `~~8~~ ~~9~~ ~~10~~ 11 ~~12~~

` 13 `~~14~~ ~~15~~ ~~16~~ 17 ~~18~~

` 19 `~~20~~ ~~21~~ ~~22~~ 23 ~~24~~

Although it can't be seen, number four is also striked through. We've identified 2 and 3 (the original primes) and eliminated 4 and 6 and then in every n-th set 2+n*6, 3+n*6, 4+n*6, 6+n*6. We then use regular sieve for the remaining numbers.

For the first three primes the code goes like this:

`fn wheelSieve nr =`

`(`

`local count = nr/30`

`local wheel = #{1,7,11,13,17,19,23,29}`

`local sieve = #{7,11,13,17,19,23,29,nr}`

`for i = 1 to count do`

`for w in wheel do`

`sieve[30*i + w] = true`

`local rootNr = nr^0.5`

`for p in sieve while p < rootNr do`

`(`

`local pSquared = p^2, twoP = 2*p`

`for i = pSquared to nr by twoP do`

`sieve[i] = false`

`)`

`sieve + #{2..3,5} - #{nr..sieve.count}`

`)`

First we define how many times the task has to be repeated – we use integer division as we will treat the first 30 differently. Then the wheel is defined (as there are just a few elements I didn't bother writing the code to calculate it for me; if using four or more primes, it would make sense, though) and then the base for our sieve which differs from the wheel only in that it doesn't include one and that it includes the upper limit itself. The reason it doesn't include the primes we were testing for in the beginning is that we would be repeating the work there, and inclusion of the upper boundary is there to preallocate the memory for the `bitarray`

.

Then it's just repeated addition of the values and sieving the potential primes. We get rid of surplus elements of the `bitarray`

, i.e. the upper value set in the beginning and the primes that might be bigger than that in the last line, together with adding the primes we've initially excluded. All in all, it's only slightly faster than `sieveEratosthenes`

but at least it isn't slower. The solution to the problem then is quite an elementary one, just add up all the primes and return the result:

`(`

`fn wheelSieve nr =`

`(`

`local count = nr/30`

`local wheel = #{1,7,11,13,17,19,23,29}`

`local sieve = #{7,11,13,17,19,23,29,nr}`

`for i = 1 to count do`

`for w in wheel do`

`sieve[30*i + w] = true`

`local rootNr = nr^0.5`

`for p in sieve while p < rootNr do`

`(`

`local pSquared = p^2, twoP = 2*p`

`for i = pSquared to nr by twoP do`

`sieve[i] = false`

`)`

`sieve + #{2..3,5} - #{nr..sieve.count}`

`)`

`local primes = wheelSieve 2000000`

`local result = 0L`

`for p in primes do`

`result += p`

`result`

`)`

DISCLAIMER: All scripts and snippets are provided as is under Creative Commons Zero (public domain, no restrictions) license. The author and this blog cannot be held liable for any loss caused as a result of inaccuracy or error within these web pages. Use at your own risk.

## This Post needs Your Comment!

Leave a Comment