Finding the n largest primes under m

Question

The algorithm below is the Sieve of Eratosthenes which I have implemented in Python. This algorithm finds the first primes up to a given value, namely n.

The way this algorithm works is by printing all primes starting from 2 up to n. What I want to do is the other way around, that is, starting from n and going down to m, or some random slice e.g the largest 100 numbers. So for example.. if i had n being some huge number (e.g. 200 million), I would be able to print the largest 100 primes from 200 million (NOT FROM THE 200TH MILLION PRIME NUMBER, BUT A PRIME THAT IS EQUAL TO/LESS THAN 200 MILLION!).

Can anyone help?

def primes(n):    
    array = [i for i in range(2,n+1)] 
    p = 2 

    while p <= n:
        i = 2*p
        while i <= n:
            array[i-2] = 0
            i += p
        p += 1  

    return [num for num in array if num > 0]

Answer 1

Here's one way:

def largest_primes_under(number, cap):
    n = cap - 1
    while number and n >= 2:
        if all(n % d for d in range(2, int(n ** 0.5 + 1))):
            yield n
            number -= 1
        n -= 1

Demo:

for p in largest_primes_under(10, 10**9):
    print(p)

Output:

999999937
999999929
999999893
999999883
999999797
999999761
999999757
999999751
999999739
999999733

Answer 2

Here's a modified version of the range sieve program I linked to in the comments. It uses the prime number theorem to get an approximation of the density of primes near hi and uses that density to estimate lo such that the number of primes in range(lo, hi) >= num. For further info please see the Mathworld article on the Prime Counting Function.

Generally, there's no simple way to accurately estimate the number of primes in a given range. FWIW, in special cases it's easy to show that there are no primes in a given range, eg every number in range(n! + 2, n! + n + 1) is clearly divisible by a number in range(2, n + 1).

The code below (hopefully :) ) over-estimates the required range. But we don't want it to waste time by using a value for lo that's too low. In some cases, lo won't be small enough, especially when num is very small, but it should be ok if you use sensible values for num.

#! /usr/bin/env python

''' Prime range sieve.

    Written by PM 2Ring 2014.10.15

    2015.05.24 Modified to find the `num` highest primes < `hi` 
'''

import sys
from math import log

def potential_primes():
    ''' Make a generator for 2, 3, 5, & thence all numbers coprime to 30 '''
    s = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
    for i in s:
        yield i
    s = (1,) + s[3:] 
    j = 30
    while True:
        for i in s:
            yield j + i
        j += 30


def range_sieve(lo, hi):
    ''' Create a list of all primes in the range(lo, hi) '''

    #Mark all numbers as prime
    primes = [True] * (hi - lo)

    #Eliminate 0 and 1, if necessary
    for i in range(lo, min(2, hi)):
        primes[i - lo] = False

    ihi = int(hi ** 0.5)
    for i in potential_primes():
        if i > ihi: 
            break

        #Find first multiple of i: i >= i*i and i >= lo
        ilo = max(i, 1 + (lo - 1) // i ) * i

        #Determine how many multiples of i >= ilo are in range
        n = 1 + (hi - ilo - 1) // i

        #Mark them as composite
        primes[ilo - lo : : i] = n * [False]

    return [i for i,v in enumerate(primes, lo) if v]


def main():
    hi = int(sys.argv[1]) if len(sys.argv) > 1 else 1000000
    num = int(sys.argv[2]) if len(sys.argv) > 2 else 1000

    #Estimate required range using the prime number theorem
    #The constant 1.25506 is from the MathWorld article on
    #the Prime Counting Function
    d = num * log(hi) * 1.25506

    #An alternate estimator. Sometimes estimates a little too low. 
    #d = num * (log(hi) + log(log(hi)))

    lo = max(2, hi - int(d))
    print 'lo =', lo

    primes = range_sieve(lo, hi)
    print len(primes), 'primes found'

    for i, p in enumerate(reversed(primes[-num:]), 1): print '%2d: %d' % (i, p)


if __name__ == '__main__':
    main()

output when called with $ python range_sieve.py 1234567890 60

lo = 1234566314
67 primes found
 1: 1234567811
 2: 1234567801
 3: 1234567799
 4: 1234567783
 5: 1234567759
 6: 1234567723
 7: 1234567717
 8: 1234567679
 9: 1234567669
10: 1234567583
11: 1234567571
12: 1234567553
13: 1234567517
14: 1234567469
15: 1234567361
16: 1234567337
17: 1234567309
18: 1234567303
19: 1234567291
20: 1234567289
21: 1234567267
22: 1234567261
23: 1234567249
24: 1234567223
25: 1234567181
26: 1234567171
27: 1234567151
28: 1234567129
29: 1234567097
30: 1234567093
31: 1234567069
32: 1234567063
33: 1234567007
34: 1234566997
35: 1234566959
36: 1234566953
37: 1234566947
38: 1234566899
39: 1234566887
40: 1234566863
41: 1234566769
42: 1234566761
43: 1234566737
44: 1234566713
45: 1234566691
46: 1234566677
47: 1234566643
48: 1234566601
49: 1234566581
50: 1234566559
51: 1234566547
52: 1234566533
53: 1234566527
54: 1234566521
55: 1234566499
56: 1234566497
57: 1234566491
58: 1234566481
59: 1234566479
60: 1234566467

Answer 3

Didn't run the code. It should work, but change where it is needed. I just posted it to give you the intuition.

def primes(n):    
    array = [1 for i in range(0,n)]
    array[0] = 0
    array[1] = 0 
    p = 2 
    while p <= n:
        i = 2*p
        while i <= n:
            array[i] = 0
            i += p
        p += 1
        while (array[p] == 0):
            p += 1 
    result = [] 
    i = 0
    while i <= n:
        if array[i]==1:
            result.append(i)
        i += 1
    return result

Answer 4

def primes(n):
    p = n/2
    return (n%i for i in range(2, p+1))

def test_primes(n, n_largest):
    ret = (k for k in range(n, 0, -1) if all(primes(k)))
    return [ret.next() for v in range(n_largest)]

usage

>>print test_primes(20, 2)
[17, 19]

to get the n largest number under

the function "primes" get the rest of the division of an element by all the element less than its middle part then "test_primes" test if all thoses rests are 0 or not. if all the rests are not 0 then this element is a prime number so we print it