Does a prime h^2+1 divide (2p)^2+1 for all prime p?

test.py

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")