Fastest way to find prime in intervals

Question

I tried searching StackOverflow and some other sources to find the fastest way to get a prime number within some interval but I didn't find any efficient way, so here is my code:

def prime(lower,upper):
   prime_num = []
   for num in range(lower, upper + 1):
       # all prime numbers are greater than 1
       if num > 1:
           for i in range(2, num):
               if (num % i) == 0:
                   break
           else:
               prime_num.append(num)
   return prime_num

Can I make this more efficient?

I tried finding my answer in Fastest way to find prime number but I didn't find prime numbers in intervals.

Answer 1

Yes. In general, use:

• Sieve of Eratosthenes
• Miller-Rabin Primality Test

Other options include trying a compiled language, like C++, but I think that isn't what you're looking for, since you've asked about Python.

Answer 2

I wrote an equation based prime finder using MillerRabin, it's not as fast as a next_prime finder that sieves, but it can create large primes and you get the equation it used to do it. Here is an example. Following that is the code:

In [5]: random_powers_of_2_prime_finder(1700)                                                                                                                               
Out[5]: 'pow_mod_p2(27667926810353357467837030512965232809390030031226210665153053230366733641224969190749433786036367429621811172950201894317760707656743515868441833458231399831181835090133016121983538940390210139495308488162621251038899539040754356082290519897317296011451440743372490592978807226034368488897495284627000283052473128881567140583900869955672587100845212926471955871127908735971483320243645947895142869961737653915035227117609878654364103786076604155505752302208115738401922695154233285466309546195881192879100630465, 2**1700-1, 2**1700) = 39813813626508820802866840332930483915032503127272745949035409006826553224524022617054655998698265075307606470641844262425466307284799062400415723706121978318083341480570093257346685775812379517688088750320304825524129104843315625728552273405257012724890746036676690819264523213918417013254746343166475026521678315406258681897811019418959153156539529686266438553210337341886173951710073382062000738529177807356144889399957163774682298839265163964939160419147731528735814055956971057054406988006642001090729179713'

or use to get the answer directly:

In [6]: random_powers_of_2_prime_finder(1700, withstats=False)                                                                                                              
Out[6]: 4294700745548823167814331026002277003506280507463037204789057278997393231742311262730598677178338843033513290622514923311878829768955491790776416394211091580729947152858233850115018443160652214481910152534141980349815095067950295723412327595876094583434338271661005996561619688026571936782640346943257209115949079332605276629723961466102207851395372367417030036395877110498443231648303290010952093560918409759519145163112934517372716658602133001390012193450373443470282242835941058763834226551786349290424923951

The code:

import random
import math

def primes_sieve2(limit):
    a = [True] * limit
    a[0] = a[1] = False

    for (i, isprime) in enumerate(a):
        if isprime:
            yield i
            for n in range(i*i, limit, i):
                a[n] = False

def ltrailing(N):
    return len(str(bin(N))) - len(str(bin(N)).rstrip('0'))


def pow_mod_p2(x, y, z):
    "4-5 times faster than pow for powers of 2"
    number = 1
    while y:
        if y & 1:
            number = modular_powerxz(number * x, z)
        y >>= 1
        x = modular_powerxz(x * x, z)
    return number

def modular_powerxz(num, z, bitlength=1, offset=0):
   xpowers = 1<<(z.bit_length()-bitlength)
   if ((num+1) & (xpowers-1)) == 0:

      return ( num & ( xpowers -bitlength)) + 2
   elif offset == -1:
      return ( num & ( xpowers -bitlength)) + 1
   elif offset == 0:
      return ( num & ( xpowers -bitlength))
   elif offset == 1:
      return ( num & ( xpowers -bitlength)) - 1
   elif offset == 2:
      return ( num & ( xpowers -bitlength)) - 2

def MillerRabin(N, primetest, iterx, powx, withstats=False): 
  primetest = pow(primetest, powx, N) 
  if withstats == True:
     print("first: ",primetest) 
  if primetest == 1 or primetest == N - 1: 
    return True 
  else: 
    for x in range(0, iterx-1): 
       primetest = pow(primetest, 2, N) 
       if withstats == True:
          print("else: ", primetest) 
       if primetest == N - 1: return True 
       if primetest == 1: return False 
  return False 

PRIMES=list(primes_sieve2(1000000))


def mr_isprime(N, withstats=False):
    if N == 2:
      return True
    if N % 2 == 0:
      return False
    if N < 2:
        return False
    if N in PRIMES:
          return True
    for xx in PRIMES:
       if N % xx == 0:
          return False
    iterx = ltrailing(N - 1)
    k = pow_mod_p2(N, (1<<N.bit_length())-1, 1<<N.bit_length()) - 1
    t = N >> iterx
    tests = [k+1, k+2, k, k-2, k-1]
    for primetest in tests:
        if primetest >= N:
            primetest %= N
        if primetest >= 2:
            if MillerRabin(N, primetest, iterx, t, withstats) == False:
                return False
    return True

def lars_last_modulus_powers_of_two(hm):
   return math.gcd(hm, 1<<hm.bit_length())


def random_powers_of_2_prime_finder(powersnumber, primeanswer=False, withstats=True):
    while True:
       randnum = random.randrange((1<<(powersnumber-1))-1, (1<<powersnumber)-1,2)
       while lars_last_modulus_powers_of_two(randnum) == 2 and  mr_isprime(randnum//2) == False:
         randnum = random.randrange((1<<(powersnumber-1))-1, (1<<powersnumber)-1,2)
       answer = randnum//2
       # This option makes the finding of a prime much longer, i would suggest not using it as 
       # the whole point is a prime answer. 
       if primeanswer == True:
          if mr_isprime(answer) == False:
            continue
       powers2find = pow_mod_p2(answer, (1<<powersnumber)-1, 1<<powersnumber)
       if mr_isprime(powers2find) == True:
          break
       else:  
          continue
    if withstats == False:
      return powers2find
    elif withstats == True:
      return f"pow_mod_p2({answer}, 2**{powersnumber}-1, 2**{powersnumber}) = {powers2find}"
    return powers2find

def nextprime(N):
   N+=2
   while not mr_isprime(N):
      N+=2
   return N

def get_primes(lower, upper):
    lower = lower|1
    upper = upper|1
    vv = []
    if mr_isprime(lower):
      vv.append(lower)
    else:
      vv=[nextprime(lower)]
    while vv[-1] < upper:
       vv.append(nextprime(vv[-1]))
    return vv

Here is an example like yours:

In [1538]: cc = get_primes(1009732533765201, 1009732533767201)                                                                                                

In [1539]: print(cc)                                                                                                                                          
[1009732533765251, 1009732533765281, 1009732533765289, 1009732533765301, 1009732533765341, 1009732533765379, 1009732533765481, 1009732533765493, 1009732533765509, 1009732533765521, 1009732533765539, 1009732533765547, 1009732533765559, 1009732533765589, 1009732533765623, 1009732533765749, 1009732533765751, 1009732533765757, 1009732533765773, 1009732533765821, 1009732533765859, 1009732533765889, 1009732533765899, 1009732533765929, 1009732533765947, 1009732533766063, 1009732533766069, 1009732533766079, 1009732533766093, 1009732533766109, 1009732533766189, 1009732533766211, 1009732533766249, 1009732533766283, 1009732533766337, 1009732533766343, 1009732533766421, 1009732533766427, 1009732533766457, 1009732533766531, 1009732533766631, 1009732533766643, 1009732533766667, 1009732533766703, 1009732533766727, 1009732533766751, 1009732533766763, 1009732533766807, 1009732533766811, 1009732533766829, 1009732533766843, 1009732533766877, 1009732533766909, 1009732533766933, 1009732533766937, 1009732533766973, 1009732533767029, 1009732533767039, 1009732533767093, 1009732533767101, 1009732533767147, 1009732533767159, 1009732533767161, 1009732533767183, 1009732533767197, 1009732533767233]

Answer 3

Yes, you can make it more efficient. As has been mentioned, for very large numbers use Miller-Rabin. For smaller ranges use the Sieve of Eratosthenes. However, apart from that, your prime checker code is very inefficient.

• 2 is the only even prime number, which can save you doing half the work you do.

• You only need to check up to the square root of the number you are testing. In any pair of factors: f and n/f, one is guaranteed to be less then or equal to the square root of the number being tested. Once you find a factor then the number is composite.

My Python is not good, so this is in pseudocode:

isPrime(num)

  // Negatives, 0, 1 are not prime.
  if (num < 2) return false

  // Even numbers: 2 is the only even prime.
  if (num % 2 == 0) return (num == 2)
  
  // Odd numbers have only odd factors.
  limit <- 1 + sqrt(num)
  for (i <- 3 to limit step 2)
    if (num % i == 0) return false
  
  // No factors found so num is prime
  return true
  
end isPrime