Search code examples
pythonlogicprimesmemory-efficient

Fastest way to find prime in intervals


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.


Solution

  • 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]