#!/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)