## Two Sieving Problems

### August 27, 2013

We have today two exercises devoted to enumerating prime numbers:

1) How many prime numbers are there between 1050 and 1050 + 106?

2) How many prime numbers between 1,000,000 and 2,000,000 are congruent to 1 (mod 4)?

Your task is to write programs that answer the two questions given above; you might note the hint in the title of the exercise. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

About these ads

Pages: 1 2 3

### 2 Responses to “Two Sieving Problems”

1. Graham said

My first answers are in J; it’s definitely cheating when your language has a
“prime factors of” function.
Problem 1:

```+/@(1=#@q:)(10x^50)+1+i.1e6%2
```

This takes a very long time to get anywhere; the second is much faster:

```+/@(1=#@q:)1e6+1+4*1e6%4
```

I originally had simpler answers, but decided that only working over odd
numbers (first question) or numbers equivalent to 1 mod 4 (second question) was
worth the uglification of my code.

My second answers are in Ruby; I can’t imagine trying to do these in a language
like C or C++, even with the help of multiprecision libraries like GMP or NTL.
I went with the Miller-Rabin primality test.

```require 'mathn'

def decomp n
s, d = 0, n
while d.even?
s, d = s + 1, d >> 1
end
return s, d
end

def modpow n, e, m
r = 1
while e > 0
if e.odd?
r = (r * n) % m
end
e >>= 1
n = (n * n) % m
end
return r
end

def miller_rabin n, k=42
s, d= decomp(n - 1)
k.times do
a = 2 + rand(n - 4)
x = modpow a, d, n
next if [1, n - 1].include? x
flag = (s - 1).times do
x = (x * x) % n
return false if x == 1
break n - 1 if x == n - 1
end
next if flag == n - 1
return false
end
return true
end

class Integer
def prime?
return false if self < 2
return true if self == 2
return false if self.even?
return miller_rabin self
end
end

def problem1
(10**50 + 1 .. 10**50 + 10**6).step(2).select(&:prime?).length
end

def problem2
(10**6 + 1 .. 2 * 10**6).step(4).select(&:prime?).length
end
```
2. Paul said

A Python version with segments.

```import time
from math import ceil, sqrt

from gmpy2 import is_prime

from ma.primegmp import sieve

def erase(flags, primes, q, n):
for pk, qk in zip(primes, q):
for i in xrange(qk, n, pk):
flags[i] = 0

def gen_sieve_seg(l, r, number_miller_rabin=25):
"""segmented sieve
uses prime numbes <= 1000000
checks candidates with miller rabin
"""
primes = list(sieve(1000000))
primes.remove(2)
q = [(-(l + 1 + pk) // 2) % pk for pk in primes]
n = ((r - l + 1) // 2)
flags = [1] * n
erase(flags, primes,q, n)
for i, f in enumerate(flags):
if f:
cand = l + i * 2 + 1
if is_prime(cand, number_miller_rabin):
yield cand

def gen_sieve_seg_4np1(l, r):
"""segmented sieve
ony numbers 4*n+1 are sieved
"""
primes = list(sieve(int(ceil(sqrt(r)))))
primes.remove(2)
assert not l % 4
q = [(-(l + 1 - (pk % 4) * pk) // 4) % pk for pk in primes]
n = (r - l) // 4
flags = [1] * n
erase(flags, primes,q, n)
for i, f in enumerate(flags):
if f:
yield l + i * 4 + 1

low = 10 ** 50
high = low + 10 ** 6

t0 = time.clock()
print sum(1 for p in gen_sieve_seg(low, high, 5)),
print "{:6.1f} sec".format(time.clock() - t0)

t0 = time.clock()
print sum(1 for p in gen_sieve_seg_4np1(1000000, 2000000)),
print "{:6.3f} sec".format(time.clock() - t0)

"""
8737    4.3 sec
35241  0.044 sec
"""
```