project-euler/518_prime_triples.py

39 lines
779 B
Python

#!/usr/bin/env python3
def getprimes(n):
sieve = [True] * n
for p in range(3, n, 2):
if not sieve[p]:
continue
for i in range(p**2, n, 2*p):
sieve[i] = False
return sieve
def is_prime(n):
if n == 2:
return True
return n % 2 and sieve[n]
N = 10**8
sieve = getprimes(N)
primes = [2] + [n for n in range(3, N, 2) if sieve[n]]
for n in range(2, int(N**0.5)+1):
b = n**2
if n % 2:
b *= 2
for i in range(b, N, b):
sieve[i] = n
s = 0
for p in primes:
d = sieve[p+1]
k = d + 1
while ((p+1) * k**2)/d**2 <= N:
if is_prime(((p+1)*k)//d - 1) and is_prime(((p+1)*k**2)//d**2 - 1):
s += p + ((p+1)*k)//d - 1 + ((p+1)*k**2)//d**2 - 1
k += 1
print(s)