p263

p263.py

from prime_sieve.array import PrimeArraySieve

def prime_factors_with_exponent(n):
    factors = []
    d = 2
    while n > 1:
        count = 0
        while n % d == 0:
            count += 1
            n /= d
        if count > 0:
            factors.append([d, count])
        d = d + 1
        if d * d > n:
            if n > 1:
                factors.append([int(n), 1])
            break
    return factors

def is_practical(x):
    pf = prime_factors_with_exponent(x)
    if pf[0][0] != 2: #Start by checking is 2 is a prime factor for completeness
        return False
    if len(pf) == 1: #This means 2 is the only prime factor so x = 2^e => it is practical 
        return True
    for x in range(1, len(pf)):
        p1, e1 = pf[x]
        total = 1
        for y in range(x):
            p, e = pf[y]
            if p != p1:
                total *= ((pow(p, (e + 1)) - 1)//(p - 1))
                #This is calculating the sum of divisors of x without the prime factor, p1
        if p1 > total + 1:
            #If a single prime fails the condition, then x is not practical
            return False  
    return True

def compute():
    limit = 10**8
    sieve = PrimeArraySieve()
    primes = sieve.primes_in_range(1, limit) #Super fast prime sieve
    total = 0
    count = 0
    x = 1
    while True:
        p1, p2, p3, p4= primes[x], primes[x+1], primes[x+2], primes[x+3]
        x += 1
        if limit - p1 < 1000:
            #If we reach the end of our primes
            #Generate a new set of 10^8 primes
            limit += 10**8
            primes = sieve.primes_in_range(p1 + 1, limit)
            x = 0
        if (p4 - p3 == 6) and (p3 - p2 == 6) and (p2 - p1 == 6):
            #Three consecutive sexy primes have been found
            n = p1 + 9
            if all([is_practical(x) for x in [n-8, n-4, n, n+4, n+8]]):
                #engineers' paradise has been found
                count += 1
                total += n
                if count == 4:
                    return int(total)