Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Does a prime h^2+1 divide (2p)^2+1 for all prime p?
from math import sqrt
max = 800000
smallest_prime_factor = []
for num in range(max): smallest_prime_factor.append(0)
for num in range(2, max):
if smallest_prime_factor[num] == 0:
smallest_prime_factor[num] = num
for num2 in range(max//num):
if smallest_prime_factor[num*num2] == 0:
smallest_prime_factor[num*num2] = num
def prime_factorization(num):
primes = {}
while num != 1:
if smallest_prime_factor[num] not in primes:
primes[smallest_prime_factor[num]] = 0
primes[smallest_prime_factor[num]] += 1
num //= smallest_prime_factor[num]
return primes
for num in range(2, 400):
if smallest_prime_factor[num] == num:
test = 2*num*2*num+1
pf_test = prime_factorization(test)
prime_found = None
for prime in pf_test:
if prime == test:
prime_found = test
break
radicand = int(sqrt(prime-1))
if radicand*radicand+1 == prime:
prime_found = prime
break
if prime_found != None: print(num, pf_test, prime_found)
else: print(num, pf_test, "FAIL")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.