No of Pairs of consecutive prime numbers having difference of 6 like (23,29) from 1 to 2 billion

Question

How to find number of pairs of consecutive prime numbers having difference of 6 like (23,29) from 1 to 2 billion (using any programming language and without using any external libraries) with considering time complexity?

1. Tried sieve of eratosthenes but getting consecutive primes is challenge

2. Used generators but time complexity is very high

The code is:

def gen_numbers(n):
    for ele in range(1,n+1):
        for i in range(2,ele//2):
            if ele%i==0:
                break
        else:
            yield ele
    prev=0
    count=0
    for i in gen_numbers(2000000000):
        if i-prev==6:
            count+=1
        prev = i

Answer 1

Interesting question! I have recently been working on Sieve of Eratosthenes prime generators. @Hans Olsson says

You should use segmented sieve to avoid memory issue: en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Segmented_sieve

I agree, and happen to have one which I hacked to solve this question. Apologies in advance for the length and non Pythonic-ness. Sample output:

$ ./primes_diff6.py 100
7 prime pairs found with a difference of 6.
( 23 , 29 ) ( 31 , 37 ) ( 47 , 53 ) ( 53 , 59 ) ( 61 , 67 ) ( 73 , 79 ) ( 83 , 89 )
25 primes found.
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 
83, 89, 97]
$ ./primes_diff6.py 1e5
1940 prime pairs found with a difference of 6.
9592 primes found.

The code:

#!/usr/bin/python -Wall
# program to find all primes smaller than n, using segmented sieve
# see https://github.com/kimwalisch/primesieve/wiki/Segmented-sieve-of-Eratosthenes

import sys

def segmentedSieve(limit):

    sqrt = int(limit ** 0.5)
    segment_size = sqrt

    prev = 0
    count = 0

    # we sieve primes >= 3
    i = 3
    n = 3

    sieve       = []
    is_prime    = [True] * (sqrt + 1)
    primes      = []
    multiples   = []
    out_primes  = []
    diff6       = []

    for low in xrange(0, limit+1, segment_size):
        sieve = [True] * segment_size

        # current segment = [low, high]
        high = min(low + segment_size -1, limit)

        # add sieving primes needed for the current segment
        #   using a simple sieve of Eratosthenese, starting where we left off
        while i * i <= high:
            if is_prime[i]:
                primes.append(i)
                multiples.append(i * i - low)

                two_i = i + i
                for j in xrange(i * i, sqrt, two_i):
                    is_prime[j] = False
            i += 2

        # sieve the current segment
        for x in xrange(len(primes)):
            k = primes[x] * 2
            j = multiples[x]

            while j < segment_size:     # NB: "for j in range()" doesn't work here.
                sieve[j] = False
                j += k

            multiples[x] = j - segment_size

        # collect results from this segment
        while n <= high:
            if sieve[n - low]:
                out_primes.append(n)
                if n - 6 == prev:
                    count += 1
                    diff6.append(n)
                prev = n
            n += 2

    print count, "prime pairs found with a difference of 6."
    if limit < 1000:
        for x in diff6:
            print "(", x-6, ",", x, ")",
        print
    return out_primes

# Driver Code
if len(sys.argv) < 2:
    n = 500
else:
    n = int(float(sys.argv[1]))

primes = [2] + segmentedSieve(n)

print len(primes), "primes found."
if n < 1000:
    print primes

This might work as-is if you run it for size 2e9 (2 billion) and subtract the result of size 1e9 (1 billion).

EDIT

Performance info, requested by @ValentinB.

$ time ./primes_diff6.py 2e9 
11407651 prime pairs found with a difference of 6. 
98222287 primes found. 

real 3m1.089s 
user 2m56.328s 
sys  0m4.656s

... on my newish laptop, 1.6 GHz i5-8265U, 8G RAM, Ubuntu on WSL, Win10

I found a mod 30 prime wheel here in a comment by Willy Good that is about 3x faster than this code at 1e9, about 2.2x faster at 2e9. Not segmented, the guts is a Python generator. I'm wondering if it could be segmented or changed to use a bit array to help its memory footprint without otherwise destroying its performance.

END EDIT

Answer 2

This will require storing all primes from 0 to sqrt(2000000000) so memory wise it's not optimal but maybe this will do for you ? You're going to have to go for a more complex sieve if you want to be more efficient though.

n = 2000000000
last_prime = 3
prime_numbers_to_test = [2, 3]
result = 0

for i in range(5, n, 2):
    for prime in prime_numbers_to_test:
        # Not prime -> next
        if i % prime == 0:
            break

    else:
        # Prime, test our condition
        if i - last_prime == 6:
            result += 1

        last_prime = i

        if i**2 < n:
            prime_numbers_to_test.append(i)

print(result)

EDIT This code yielded a result of 11,407,651 pairs of consecutive primes with a difference of 6 for n = 2,000,000,000

Answer 3

Here's a demonstration along the lines of what I understood as user448810's intention in their comments. We use the primes to mark off, meaning sieve, only relevant numbers in the range. Those are numbers of the form 6k + 1 and 6k - 1.

Python 2.7 code:

# https://rosettacode.org/wiki/Modular_inverse
def extended_gcd(aa, bb):
    lastremainder, remainder = abs(aa), abs(bb)
    x, lastx, y, lasty = 0, 1, 1, 0
    while remainder:
        lastremainder, (quotient, remainder) = remainder, divmod(lastremainder, remainder)
        x, lastx = lastx - quotient*x, x
        y, lasty = lasty - quotient*y, y
    return lastremainder, lastx * (-1 if aa < 0 else 1), lasty * (-1 if bb < 0 else 1)

def modinv(a, m):
    g, x, y = extended_gcd(a, m)
    if g != 1:
        raise ValueError
    return x % m


from math import sqrt, ceil

n = 2000000000
sqrt_n = int(sqrt(n))

A = [True] * (sqrt_n + 1)

for i in xrange(2, sqrt_n // 2):
  A[2*i] = False

primes = [2]

for i in xrange(3, sqrt_n, 2):
  if A[i]:
    primes.append(i)
    c = i * i
    while c <= sqrt_n:
      A[c] = False
      c = c + i 

print "Primes with a Difference of 6"
print "\n%s initial primes to work from." % len(primes)

lower_bound = 1000000000
upper_bound = 1000001000
range = upper_bound - lower_bound

print "Range: %s to %s" % (lower_bound, upper_bound)

# Primes of the form 6k + 1

A = [True] * (range // 6 + 1)

least_6k_plus_1 = int(ceil((lower_bound - 1) / float(6)))
most_6k_plus_1 = (upper_bound - 1) // 6

for p in primes[2:]:
  least = modinv(-6, p)
  least_pth_over = int(least + p * ceil((least_6k_plus_1 - least) / float(p)))
  c = int(least_pth_over - least_6k_plus_1)
  while c < len(A):
    A[c] = False
    c = c + p

print "\nPrimes of the form 6k + 1:"

for i in xrange(1, len(A)):
  if A[i] and A[i - 1]:
    p1 = (i - 1 + least_6k_plus_1) * 6 + 1
    p2 = (i + least_6k_plus_1) * 6 + 1
    print p1, p2


# Primes of the form 6k - 1

A = [True] * (range // 6 + 1)

least_6k_minus_1 = int(ceil((lower_bound + 1) / float(6)))
most_6k_minus_1 = (upper_bound + 1) // 6

for p in primes[2:]:
  least = modinv(6, p)
  least_pth_over = int(least + p * ceil((least_6k_minus_1 - least) / float(p)))
  c = int(least_pth_over - least_6k_minus_1)
  while c < len(A):
    A[c] = False
    c = c + p

print "\nPrimes of the form 6k - 1:"

for i in xrange(1, len(A)):
  if A[i] and A[i - 1]:
    p1 = (i - 1 + least_6k_minus_1) * 6 - 1
    p2 = (i + least_6k_minus_1) * 6 - 1
    print p1, p2

Output:

Primes with a Difference of 6

4648 initial primes to work from.
Range: 1000000000 to 1000001000

Primes of the form 6k + 1:
1000000087 1000000093
1000000447 1000000453
1000000453 1000000459

Primes of the form 6k - 1:
1000000097 1000000103
1000000403 1000000409
1000000427 1000000433
1000000433 1000000439
1000000607 1000000613

In order to count consecutive primes, we have to take into account the interleaving lists of primes 6k + 1 and 6k - 1. Here's the count:

# https://rosettacode.org/wiki/Modular_inverse
def extended_gcd(aa, bb):
    lastremainder, remainder = abs(aa), abs(bb)
    x, lastx, y, lasty = 0, 1, 1, 0
    while remainder:
        lastremainder, (quotient, remainder) = remainder, divmod(lastremainder, remainder)
        x, lastx = lastx - quotient*x, x
        y, lasty = lasty - quotient*y, y
    return lastremainder, lastx * (-1 if aa < 0 else 1), lasty * (-1 if bb < 0 else 1)

def modinv(a, m):
    g, x, y = extended_gcd(a, m)
    if g != 1:
        raise ValueError
    return x % m


from math import sqrt, ceil
import time

start = time.time()

n = 2000000000
sqrt_n = int(sqrt(n))

A = [True] * (sqrt_n + 1)

for i in xrange(2, sqrt_n // 2):
  A[2*i] = False

primes = [2]

for i in xrange(3, sqrt_n, 2):
  if A[i]:
    primes.append(i)
    c = i * i
    while c <= sqrt_n:
      A[c] = False
      c = c + i 

lower_bound = 1000000000
upper_bound = 2000000000
range = upper_bound - lower_bound

A = [True] * (range // 6 + 1)

least_6k_plus_1 = int(ceil((lower_bound - 1) / float(6)))
most_6k_plus_1 = (upper_bound - 1) // 6

for p in primes[2:]:
  least = modinv(-6, p)
  least_pth_over = int(least + p * ceil((least_6k_plus_1 - least) / float(p)))
  c = int(least_pth_over - least_6k_plus_1)
  while c < len(A):
    A[c] = False
    c = c + p

B = [True] * (range // 6 + 1)

least_6k_minus_1 = int(ceil((lower_bound + 1) / float(6)))
most_6k_minus_1 = (upper_bound + 1) // 6

for p in primes[2:]:
  least = modinv(6, p)
  least_pth_over = int(least + p * ceil((least_6k_minus_1 - least) / float(p)))
  c = int(least_pth_over - least_6k_minus_1)
  while c < len(B):
    B[c] = False
    c = c + p

total = 0

for i in xrange(1, max(len(A), len(B))):
  if A[i] and A[i - 1] and not B[i]:
    total = total + 1
  if B[i] and B[i - 1] and not A[i - 1]:
    total = total + 1

# 47374753 primes in range 1,000,000,000 to 2,000,000,000
print "%s consecutive primes with a difference of 6 in range %s to %s." % (total, lower_bound, upper_bound)
print "--%s seconds" % (time.time() - start)

Output:

5317860 consecutive primes with a difference of 6 in range 1000000000 to 2000000000.
--193.314619064 seconds

Answer 4

Python isn't the best language to write this in, but since that's what we're all doing...

This little segmented sieve finds the answer 5317860 in 3:24

import math

# Find primes < 2000000000

sieve = [True]*(int(math.sqrt(2000000000))+1)

for i in range(3,len(sieve)):
    if (sieve[i]):
        for j in range(2*i, len(sieve), i):
            sieve[j] = False

smallOddPrimes = [i for i in range(3,len(sieve),2) if sieve[i]]

# Check primes in target segments

total=0
lastPrime=0

for base in range(1000000000, 2000000000, 10000000):
    sieve = [True]*5000000
    for p in smallOddPrimes:
        st=p-(base%p)
        if st%2==0: #first odd multiple of p
            st+=p
        for i in range(st//2,len(sieve),p):
            sieve[i]=False
    for prime in [i*2+base+1 for i in range(0,len(sieve)) if sieve[i]]:
        if prime == lastPrime+6:
            total+=1
        lastPrime = prime

print(total)