Reputation: 29
Goldbach Conjecture in Python
I have attempted to write a code that returns a single pair that satisfies the Goldbach Conjecture for a given N. The conjecture states that every even number greater than 4 can be expressed as the sum of two prime numbers. The function returns a pair that is just slightly off, for example, goldbach(34) returns (5, 31) rather than the correct answer (3, 31). Similarly, goldbach(38) returns (11, 31). Any ideas where I am going wrong here? I understand that this code is not massively efficient, however this is the way I have been asked to code for my assignment.
def eratosthenes(n):
primes = list (range(2, n+1))
for i in primes:
j=2
while i*j<= primes[-1]:
if i*j in primes:
primes.remove(i*j)
j=j+1
return primes
def odd_primes(N):
oddprimes = eratosthenes(N)
oddprimes.remove(2)
return(oddprimes)
def goldbach(N):
x, y = 0, 0
result = 0
if N % 2 == 0:
prime = odd_primes(N)
while result != N:
for i in range(len(prime)):
x = prime[i]
if result == N: break
for j in range(len(prime)):
y = prime[j]
result = x + y
if result == N: break
return x, y
Upvotes: 1
Views: 12664
Answers (9)
Reputation: 189
Here is your code fixed a bit. You already had the primes so I used them. It just checks if N minus this prime is prime.
def eratosthenes(n):
primes = list (range(2, n+1))
for i in primes:
j=2
while i*j<= primes[-1]:
if i*j in primes:
primes.remove(i*j)
j=j+1
return primes
def goldbach(N):
if N % 2 == 0 and N > 3:
prime = eratosthenes(N)
for i in prime:
if N-i in prime:
return i,N-i
return None
for i in range(4,50,2):
x,y = goldbach(i)
print(f"{i} = {x} + {y}"
else:
i = 12340
x,y = goldbach(i)
print(f"{i} = {x} + {y}")
Output:
4 = 2 + 2
6 = 3 + 3
8 = 3 + 5
10 = 3 + 7
12 = 5 + 7
14 = 3 + 11
16 = 3 + 13
18 = 5 + 13
20 = 3 + 17
22 = 3 + 19
24 = 5 + 19
26 = 3 + 23
28 = 5 + 23
30 = 7 + 23
32 = 3 + 29
34 = 3 + 31
36 = 5 + 31
38 = 7 + 31
40 = 3 + 37
42 = 5 + 37
44 = 3 + 41
46 = 3 + 43
48 = 5 + 43
12340 = 11 + 12329
Don't check for prime unless forced.
10,000,000,000,000 less than a second
This takes milliseconds for trillions. Using sympy, gmpy2 or home grown IsPrime. Just don't call for a Prime check unless you have to. For very large numbers you only call IsPrime() about 50 times.
#from sympy import isprime as IsPrime
def IsPrime(n):
if (n <= 3): return n > 1
if n%6 not in [1,5]: 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: return 2,2
elif IsPrime(number-3): return 3,number-3
else: # just 6k±1
for p,q in ((i,i+2) for i in range(5,number//2+1,6)):
if IsPrime(p) and IsPrime(number-p):
return p,number-p
if IsPrime(q) and IsPrime(number-q):
return q,number-q
return 0,0 # failed
print(goldbach(3325581707333960528))
Output:
$ python goldbach_little.py
(9781, 3325581707333950747)
Though I think the sieve solutions are all limited for large 2^100 size numbers, here is a faster implementation:
def eratosthenes(n):
sieve = [False, True] * (n // 2) # No evens
sieve[1] = False # odd non prime
sieve[2] = True # even prime
sieve[9::3] = [False] * len(sieve[9::3]) # clear threes
for p,q in ((i,i+2) for i in range(5,int(n**.5)+1,6)):
if sieve[p]:
sieve[p*p::2*p] = [False] * len(sieve[p*p::2*p])
if sieve[q]:
sieve[q*q::2*q] = [False] * len(sieve[q*q::2*q])
return sieve
def goldbach(N):
if N % 2 == 0 and N > 3:
if N == 4: return 2,2
prime = eratosthenes(N)
if prime[N-3]: return 3,N-3
for p,q in ((i,i+2) for i in range(5,N//2+1,6)):
if prime[p] and prime[N-p]: return p,N-p
if prime[q] and prime[N-q]: return q,N-q
return None
for i in range(400,550,2):
x,y = goldbach(i)
print(f"{i} = {x} + {y}")
else:
i = 123456780
x,y = goldbach(i)
print(f"{i} = {x} + {y}")
Upvotes: 0
Reputation: 151
You want to get, for any even number N, a pair of prime number whose sum is N.
For this, we will simply go through each number between 2 and N. Let's call this number i. We will then check if i and N-i are both prime.
If that is the case, then we have a pair of prime number so that their sum is N.
Indeed, i + (N-i) = N.
def is_prime(number):
for div in range(2, number//2 + 1):
if number % div == 0:
return False
return True
def compute_goldbach_pair(number):
for p1 in range(2, number):
p2 = number - p1
if is_prime(p1) and is_prime(p2):
return [p1, p2]
In this code, we return as soon as we have found a pair. But multiple pairs can exists, in this case, you'll have to modify the code a bit to collect each existing pairs.
Upvotes: 1
Reputation: 11
def is_prime(a):
for div in range(2, a//2 + 1):
if a % div == 0:
return False
return True
def goldbach(n):
summ = 0
for p1 in range(2, n):
if is_prime(p1) and is_prime(n-p1):
p2 = n-p1
return str(n)+"="+str(p1)+"+"+str(p2)
Upvotes: 1
Reputation: 4988
The Goldbach decomposition of an even number into a pair of primes numbers is not (always) unique so you need a method to select one possible solution (as required from the question). A natural possibility would be using min/max function.
I did not use the function odd_primes (just a slice) and reimplemented goldbach function using a combinatorics approach.
import itertools as it
def eratosthenes(n):
primes = list (range(2, n+1))
for p in primes:
j = 2
while p*j <= primes[-1]:
if p*j in primes:
primes.remove(p*j)
j += 1
return primes
def goldbach(n):
if n % 2 != 0 or n <= 2: raise Exception('Conjecture assumptions not satisfied.')
ps = eratosthenes(n)[1:]
# trivial cases to overcome limitations of the combinations
if n == 4: return [(2, 2)]
if n == 6: return [(3, 3)]
return max(it.combinations(ps, 2), key=lambda p: p[1] if sum(p) == n else 0)
for n in range(6, 40, 2):
print(n, goldbach(n))
Output
8 (3, 5)
10 (3, 7)
12 (5, 7)
14 (3, 11)
16 (3, 13)
18 (5, 13)
20 (3, 17)
22 (3, 19)
24 (5, 19)
26 (3, 23)
28 (5, 23)
30 (7, 23)
32 (3, 29)
34 (3, 31)
36 (5, 31)
38 (7, 31)
Upvotes: 0
Reputation: 179
It is possible to use a sieve that separates the primes 1 (mod 6) with the ones -1 (mod 6) in this way it is possible to speed up using numpy without using the for loop:
import numpy as np
def sieve(n):
P5m6 =np.ones((n//6+1), dtype=bool)
P1m6 =np.ones((n//6+1), dtype=bool)
for i in range(1,int((n**0.5+1)/6)+1):
if P5m6[i]:
P5m6[6*i*i::6*i-1]=[False]
P1m6[6*i*i-2*i::6*i-1]=[False]
if P1m6[i]:
P5m6[6*i*i::6*i+1]=[False]
P1m6[6*i*i+2*i::6*i+1]=[False]
return P5m6,P1m6
def goldbach(n):
P5m6,P1m6=sieve(n)
nmod6=n%6
if nmod6==0:
r=(6*np.nonzero(P5m6[1:n//6]&P1m6[n//6-1:0:-1])[0]+5)[0]
elif nmod6==2:
if P5m6[n//6]:
r=3
else:
r=(6*np.nonzero(P1m6[1:(n-6)//12+(n//6-1)%2+1]&P1m6[n//6-1:(n-6)//12:-1])[0]+7)[0]
elif nmod6==4:
if P1m6[n//6]:
r=3
else:
r=(6*np.nonzero(P5m6[1:n//12+(n//6)%2+1]&P5m6[n//6:n//12:-1])[0]+5)[0]
else:
r=0
return r,n-r
Upvotes: 1
Reputation: 68
Goldbach Conjecture is observed only in even numbers greater than 2, therefore I improve on Punnie's answer to catch odd numbers and any number less done 2.
def isPrime(n):
for i in range(2,n):
if n%i == 0:
return 0
return 1
def goldbach(no):
if no%2 !=0 :
return "Error {} is not an even number ".format(no)
elif no <= 2:
return "Error {} is not greater than 2, Goldbach Conjecture is observed only in even numbers greater than 2".format(no)
else:
for i in range(3,no):
if isPrime(i) == 1:
for l in range(i,no):
if isPrime(l) == 1:
if no == (i+l):
print(i,"+",l,"=",no)
no = int(input("Enter Even Number greater than 2: "))
goldbach(no) ```
Upvotes: 0
Reputation: 1
def isPrime(n):
for i in range(2,n):
if n%i == 0:
return 0
return 1
no = int(input("Enter Number: "))
for i in range(3,no):
if isPrime(i) == 1:
for l in range(i,no):
if isPrime(l) == 1:
if no == (i+l):
print(i,"+",l,"=",no)
Upvotes: 0
Reputation: 332
def eratosthenes(n):
primes = list (range(2, n+1))
for i in primes:
j=2
while i*j<= primes[-1]:
if i*j in primes:
primes.remove(i*j)
j=j+1
return primes
def odd_primes(N):
oddprimes = eratosthenes(N)
oddprimes.remove(2)
return(oddprimes)
def goldbach(N):
x, y = 0, 0
result = 0
if N % 2 == 0:
prime = odd_primes(N)
while result != N:
for i in range(len(prime)):
if result == N: break
x = prime[i]
for j in range(len(prime)):
y = prime[j]
result = x + y
if result == N: break
return x, y
is the correct version. When you have found a pair, you set x to the next prime.
Upvotes: 3
Reputation: 4343
You are assigning x before breaking the loop once your condition is met. Simply invert your break lines in the first for loop:
def goldbach(N):
x, y = 0, 0
result = 0
if N % 2 == 0:
prime = odd_primes(N)
while result != N:
for i in range(len(prime)):
if result == N: break # this line first
x = prime[i] # this line after
for j in range(len(prime)):
y = prime[j]
result = x + y
if result == N: break
return x, y
Upvotes: 4
Related Questions
- Coding Lemoine's Conjecture in Python
- Goldbach conjecture exercise (c)
- Checking Goldbach's conjecture holds up to N
- Weak Goldbach Conjecture in python
- Goldbach Python Output
- Goldbach's Conjecture - Find the number of ways an even number can be written as the sum of two primes
- Implementing Goldbach's conjecture in Haskell, lots of restrictions
- Implementing Goldbach's Conjecture in C code
- Implement Goldbach's conjecture in Haskell
- Goldbach Conjecture algorithm returning one incorrect tuple each time