Гипотеза Гольдбаха в Python
Я попытался написать код, который возвращает одну пару, которая удовлетворяет гипотезе Гольдбаха для данного N. Предположение утверждает, что каждое четное число, большее 4, может быть выражено как сумма двух простых чисел. Функция возвращает пару, которая слегка отклонена, например, goldbach(34) возвращает (5, 31), а не правильный ответ (3, 31). Аналогично, goldbach(38) возвращает (11, 31). Есть идеи, где я иду не так? Я понимаю, что этот код не очень эффективен, однако именно так меня просили написать код для моего задания.
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
9 ответов
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
это правильная версия. Когда вы нашли пару, вы устанавливаете x для следующего простого числа.
Вы назначаете x до разрыва цикла, как только ваше условие выполнено. Просто инвертируй свой break строки в первом for цикл:
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
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)
Разложение Гольдбаха четного числа на пару простых чисел не (всегда) уникально, поэтому вам нужен метод для выбора одного возможного решения (как требуется из вопроса). Естественной возможностью было бы использование
min/
max функция.
Я не пользовался функцией
odd_primes (просто кусочек) и повторно реализовал
goldbach функция с использованием комбинаторического подхода.
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))
Выход
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)
Можно использовать сито, которое отделяет простые числа 1 (mod 6) от единиц -1 (mod 6), таким образом, можно ускорить использование numpy без использования цикла for:
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
улучшение ответа Панни для обработки нечетных чисел
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)
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 Number: "))
goldbach(no)
Не проверяйте премьеру, если это не принудительно.
Для триллионов это занимает миллисекунды. Использование gmpy2 или отечественного IsPrime. Только не требуйте проверки Prime, если в этом нет необходимости. На очень большие номера вы звоните всего около 50 раз.
# 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
elif IsPrime(number - 3):
print(f"\n3 + {number-3:,} = {number}\n")
return
else:
for p in range(5, number, 6): # just odds after 3
if IsPrime(p) and IsPrime(number-p):
print(f"\n{p:,} + {number-p:,} = {N:,}\n")
return
elif IsPrime(p+2) and IsPrime(number-(p+2)):
print(f"\n{p+2:,} + {number-(p+2):,} = {N:,}\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)')
Работает очень быстро с GMPY2.is_prime():
python goldbach.py 4444444444444444
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.
107 + 4,444,444,444,444,337 = 4,444,444,444,444,444
55 function calls in 0.001 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.001 0.001 <string>:1(<module>)
1 0.000 0.000 0.001 0.001 goldbach.py:16(goldbach)
1 0.000 0.000 0.001 0.001 {built-in method builtins.exec}
1 0.001 0.001 0.001 0.001 {built-in method builtins.print}
50 0.000 0.000 0.000 0.000 {built-in method gmpy2.gmpy2.is_prime}
1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
С домашним IsPrime():
python goldbach.py 4444444444444444
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.
107 + 4,444,444,444,444,337 = 4,444,444,444,444,444
55 function calls in 3.957 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 3.957 3.957 <string>:1(<module>)
1 0.000 0.000 3.957 3.957 goldbach.py:16(goldbach)
50 3.955 0.079 3.955 0.079 goldbach.py:6(IsPrime)
1 0.000 0.000 3.957 3.957 {built-in method builtins.exec}
1 0.001 0.001 0.001 0.001 {built-in method builtins.print}
1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
10 триллионов с использованием всего Python(без GMPY2)
python goldbach.py 10000000000000
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.
29 + 9,999,999,999,971 = 10,000,000,000,000
20 function calls in 0.160 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.160 0.160 <string>:1(<module>)
1 0.000 0.000 0.160 0.160 goldbach.py:16(goldbach)
15 0.159 0.011 0.159 0.011 goldbach.py:6(IsPrime)
1 0.000 0.000 0.160 0.160 {built-in method builtins.exec}
1 0.001 0.001 0.001 0.001 {built-in method builtins.print}
1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
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)
Вы хотите получить для любого четного числа N пару простых чисел, сумма которых равна N.
Для этого мы просто пройдемся по каждому числу от 2 до N. Назовем это число i. Затем мы проверим, являются ли i и Ni простыми.
Если это так, то у нас есть пара простых чисел, сумма которых равна N.
Действительно, i + (Ni) = 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]
В этом коде мы возвращаемся, как только находим пару. Но может существовать несколько пар, в этом случае вам придется немного изменить код, чтобы собрать все существующие пары.