euler/python/e087.py

from lib_prime import primes


def euler_087():
    n_max = 50000000
    ps = primes(10**5)
    found = set()
    c = 0
    for p4 in ps:
        p4s = p4 ** 4
        if p4 > n_max:
            break
        for p3 in ps:
            p3s = p3 ** 3
            if p4s + p3s > n_max:
                break
            for p2 in ps:
                p2s = p2 ** 2
                p = p2s + p3s + p4s
                if p < n_max and not p in found:
                    c += 1
                    found.add(p)
    return c


if __name__ == "__main__":
    solution = euler_087()
    print("e087.py: " + str(solution))
    assert(solution == 1097343)