5
\$\begingroup\$

I am doing an exercise on finding the nearest pair of prime numbers p₁ and p₂ such that:

  • p₁ ≤ p₂
  • p₁ + p₂ = n (for 4 ≤ n ≤ 10⁷)
  • p₁ and p₂ are primes
  • p₂ - p₁ is minimal

I think it should work perfectly except that it got a TLE. Can someone point out how can I avoid it?

import math
import random

def _try_composite(a, d, n, s):
    if pow(a, d, n) == 1:
        return False
    for i in range(s):
        if pow(a, 2 ** i * d, n) == n - 1:
            return False
    return True  # n  is definitely composite


def is_prime(n, _precision_for_huge_n=16):
    if n in _known_primes:
        return True
    if n in (0, 1) or any((n % p) == 0 for p in _known_primes):
        return False
    d, s = n - 1, 0
    while not d % 2:
        d, s = d >> 1, s + 1
    # Returns exact according to http://primes.utm.edu/prove/prove2_3.html
    if n < 1373653:
        return not any(_try_composite(a, d, n, s) for a in (2, 3))
    if n < 25326001:
        return not any(_try_composite(a, d, n, s) for a in (2, 3, 5))


_known_primes = [2, 3]
_known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]

def main():
    n = int(input())
    if n > 4:
        for smaller in range(n // 2, -1, -1):
            if n - smaller >= smaller:
                if is_prime(n - smaller) and is_prime(smaller):
                    print (smaller, n - smaller)
                    flag = True
                    break
    else:
        print ('2 2')


main()
\$\endgroup\$
5
  • \$\begingroup\$ Since I need to minimize p_2 - p_1, it should be from n/2 to 3 and thus doesn't change? \$\endgroup\$
    – Gareth Ma
    Mar 13, 2018 at 2:29
  • \$\begingroup\$ Also using a sieve is a good idea... except sieving from 1 to 10^7 takes around 19.3 second of my code. \$\endgroup\$
    – Gareth Ma
    Mar 13, 2018 at 2:29
  • \$\begingroup\$ Let us continue this discussion in chat. \$\endgroup\$
    – Gareth Ma
    Mar 13, 2018 at 9:57
  • \$\begingroup\$ @TobySpeight Sorry, I did not copy it. I typed it out, so it is a typo. \$\endgroup\$
    – Gareth Ma
    Mar 15, 2018 at 4:43
  • 1
    \$\begingroup\$ Typing it yourself is the form of copying that I was thinking of - thanks for correcting the description. It's good to know it's not a bug in the code! \$\endgroup\$ Mar 15, 2018 at 8:42

2 Answers 2

4
\$\begingroup\$

My solution consists of a single change:

  • Changed the is_prime function into a prime sieve (from this answer) + lookup

My final code (Also handled odd numbers just for fun, definitely not because I misread the question):

def main():
    import itertools
    import math
    izip = itertools.zip_longest
    chain = itertools.chain.from_iterable
    compress = itertools.compress
    def rwh_primes2_python3(n):
        """ Input n>=6, Returns a list of primes, 2 <= p < n """
        zero = bytearray([False])
        size = n//3 + (n % 6 == 2)
        sieve = bytearray([True]) * size
        sieve[0] = False
        for i in range(int(n**0.5)//3+1):
          if sieve[i]:
            k=3*i+1|1
            start = (k*k+4*k-2*k*(i&1))//3
            sieve[(k*k)//3::2*k]=zero*((size - (k*k)//3 - 1) // (2 * k) + 1)
            sieve[  start ::2*k]=zero*((size -   start  - 1) // (2 * k) + 1)
        ans = [2,3]
        poss = chain(izip(*[range(i, n, 6) for i in (1,5)]))
        ans.extend(compress(poss, sieve))
        return ans

    string2 = "Impossible"
    n = int(input())
    sieve = [t for t in rwh_primes2_python3(n) if t <= math.floor((n // 2) / 2) * 2 + 1][::-1]
    another = [t for t in rwh_primes2_python3(n) if t >= math.floor((n // 2) / 2) * 2 + 1]

    if n > 5 and n % 2 == 0:
        for smaller in sieve:
            if n - smaller in another:
                print (smaller, n - smaller)
                break
    elif n % 2 == 1 and n != 5:
        if n - 2 in another or n-2 in sieve: print (2, n-2)
        else: print ('Impossible')
    elif n == 4:
        print ('2 2')
    else:
        print ('2 3')

main()
\$\endgroup\$
1
  • 1
    \$\begingroup\$ One advantage of this approach is that if you want to test multiple n, you just sieve for the largest one, and that cost is amortized over all the inputs. I should point out that you invoke rwh_primes2_python3(n) twice - I think you can call it just once, and then split the result into upper and lower halves (actually just one half, for the for smaller in another loop - assuming primes get sparser, that's a better choice - and we can keep sieve complete for the if n-smaller in sieve test). \$\endgroup\$ Mar 14, 2018 at 8:39
1
\$\begingroup\$

I opted to not create a sieve or a list of primes. Instead I check each prime only when forced to check. As a result, for numbers larger than 2000 many fewer primes are found. This increases the speed exponentially the larger the number becomes. 18 and 20 digit numbers require fewer than 100 prime checks usually. This code will find 18 digit even numbers pair of primes in about 20 seconds. 9 digits takes less than a second.

# Goldbach's Conjecture tester. 
from gmpy2 import is_prime
import sys
import cProfile
# or use home grown IsPrime
def IsPrime(n): 
        if (n == 2 or n == 3): 
            return True 
        if (n <= 1 or n % 2 == 0 or n % 3 == 0): 
            return False 
        for i  in range(5, int(n**.5)+1,  6): 
            if (n % i == 0 or n % (i + 2) == 0): 
                return False 
        return True 

def goldbach(number):
    if number == 4:
        print("\n2 + 2 = 4\n")
        return
    else:
        mid = int(number//2) + (1 - (number//2%2) )
        for p in range(mid, 0, -2): # just odds 
            if IsPrime(p) and IsPrime(number-p): 
                print(f"\n{p:,} + {number-p:,} = {number:,}\n")
                return
            elif IsPrime(p+2) and IsPrime(number-(p+2)):
                print(f"\n{p+2:,} + {number-(p+2):,} = {number:,}\n")
                return
        raise Exception(f"Found a counter-example to the Goldbach conjecture: {number}")

if __name__=="__main__":
    N = 1
    args = len(sys.argv)
    if args > 1:
        N = int(sys.argv[1])
    print("This is a test of Goldbach's Conjecture that for all even integers")
    print("greater than 2 there are two primes that add up to that even number.\n")
    while (N < 3 or N%2):
        N = int(input("Please enter an even number > 3 to check with Goldbach's Conjecture> "))
    cProfile.run('goldbach(N)')

Output:

python3 goldbach_near.py 
This is a test of Goldbach's Conjecture that for all even integers
greater than 2 there are two primes that add up to that even number.

Please enter an even number > 3 to check with Goldbach's Conjecture> 12345678

6,172,799 + 6,172,879 = 12,345,678

$ python3 goldbach_near.py 12348
This is a test of Goldbach's Conjecture that for all even integers
greater than 2 there are two primes that add up to that even number.


6,151 + 6,197 = 12,348
$ python3 goldbach_near.py 44000000000000000
This is a test of Goldbach's Conjecture that for all even integers
greater than 2 there are two primes that add up to that even number.


21,999,999,999,999,979 + 22,000,000,000,000,021 = 44,000,000,000,000,000

         30 function calls in 14.669 seconds
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.