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]))
```