Euler 221 – Alexandrian Numbers

When I came upon this Euler problem, and saw it was based on my namesake, I couldn’t help going for it. As with many of the Project Euler problems, the trick is to find a technique that reduces the brute force solution into something more efficient.

The functions from the original problem can be re-written to:
A = p(p + d)((p^2 + 1)/d + p)

The Alexandrian numbers are calculated from the new function, where d runs over divisors of p^2+1 and p runs over all positive integers. Used a fast Divisor routine that I Googled, then started building a Set (which is speed optimized to avoid duplicates) until length is 15K. Convert to a List, then sort and display 150,000th number. Numbers 1 thru 6 also listed to verify computation working. Runs in about 10 mins on my quad core iMac.

# e221.py

from math import sqrt
from functools import reduce


def appendEs2Sequences(sequences, es):
    result = []
    if not sequences:
        for e in es:
            result.append([e])
    else:
        for e in es:
            result += [seq + [e] for seq in sequences]
    return result


def primefactors(n):
    i = 2
    while i <= sqrt(n):
        if n % i == 0:
            l = primefactors(n // i)
            l.append(i)
            return l
        i += 1
    return [n]


def factorGenerator(n):
    p = primefactors(n)
    factors = {}
    for p1 in p:
        try:
            factors[p1] += 1
        except KeyError:
            factors[p1] = 1
    return factors


def divisors(n):
    factors = factorGenerator(n)
    divisors = []
    listexponents = [list(map(lambda x:k ** x, list(range(0, factors[k] + 1)))) for k in factors.keys()]
    listfactors = reduce(appendEs2Sequences, listexponents, [])
    for f in listfactors:
        divisors.append(reduce(lambda x, y: x * y, f, 1))
    divisors.sort()
    return divisors


if __name__ == "__main__":
    A = set()
    p = 1
    while len(A) <= 1500001:
        for d in divisors((p * p) + 1):
            A.add(int(p * (p + d) * (p + ((p * p) + 1) / d)))
        p += 1
    L = list(A)
    L.sort()
    for i in [1,2,3,4,5,6,150000]:
        print (str(i) + ": " + str(L[i-1]))