Problem 221: Alexandrian Integers

Project Euler Problem 221 Solution

EulerSolve provides an optimized solution for Project Euler Problem 221, Alexandrian Integers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary An Alexandrian integer is a positive integer obtained from an integer triple satisfying $$\frac{1}{x}+\frac{1}{y}+\frac{1}{z}=\frac{1}{xyz}.$$ For positive Alexandrian values one may choose the generating triple with exactly one negative entry, and the reported number is $A=\lvert xyz\rvert$. Problem 221 asks for the $150000$-th term when all such positive integers are listed in increasing order. The crucial point is that the search is not over arbitrary triples. The reciprocal identity collapses to a rigid divisor condition, so every valid Alexandrian integer comes from factoring a number of the form $p^2+1$. That is the mathematical reason the implementations can solve the problem without brute-forcing three independent variables. Mathematical Approach The derivation used by the implementations is a complete parameterization of the Alexandrian integers that can appear in the ordered sequence. Rewriting the reciprocal identity Take a valid triple in the form $(-p,q,r)$ with $p,q,r>0$. Substituting it into the defining identity gives $$-\frac{1}{p}+\frac{1}{q}+\frac{1}{r}=\frac{1}{(-p)qr}=-\frac{1}{pqr}.$$ Multiplying by $pqr$ removes the denominators and leaves the Diophantine equation $$qr-pq-pr=1.$$ So the reciprocal statement is equivalent to a very structured bilinear equation in the three integers....

Detailed mathematical approach

Problem Summary

An Alexandrian integer is a positive integer obtained from an integer triple satisfying

$$\frac{1}{x}+\frac{1}{y}+\frac{1}{z}=\frac{1}{xyz}.$$

For positive Alexandrian values one may choose the generating triple with exactly one negative entry, and the reported number is $A=\lvert xyz\rvert$. Problem 221 asks for the $150000$-th term when all such positive integers are listed in increasing order.

The crucial point is that the search is not over arbitrary triples. The reciprocal identity collapses to a rigid divisor condition, so every valid Alexandrian integer comes from factoring a number of the form $p^2+1$. That is the mathematical reason the implementations can solve the problem without brute-forcing three independent variables.

Mathematical Approach

The derivation used by the implementations is a complete parameterization of the Alexandrian integers that can appear in the ordered sequence.

Rewriting the reciprocal identity

Take a valid triple in the form $(-p,q,r)$ with $p,q,r>0$. Substituting it into the defining identity gives

$$-\frac{1}{p}+\frac{1}{q}+\frac{1}{r}=\frac{1}{(-p)qr}=-\frac{1}{pqr}.$$

Multiplying by $pqr$ removes the denominators and leaves the Diophantine equation

$$qr-pq-pr=1.$$

So the reciprocal statement is equivalent to a very structured bilinear equation in the three integers.

From a Diophantine equation to divisor pairs

Add $p^2$ to both sides of the previous relation:

$$qr-pq-pr+p^2=p^2+1.$$

The left-hand side factors immediately:

$$(q-p)(r-p)=p^2+1.$$

Now set

$$d=q-p,\qquad e=r-p.$$

Then $de=p^2+1$. Moreover $d$ and $e$ must be positive. If, for example, $q\le p$, then $qr-pq-pr\le pr-pq-pr=-pq\lt 0$, contradicting $qr-pq-pr=1$. Therefore every valid triple is of the form

$$(-p,\;p+d,\;p+e)\qquad\text{with}\qquad de=p^2+1.$$

The corresponding Alexandrian integer is

$$A=\lvert (-p)(p+d)(p+e)\rvert=p(p+d)(p+e).$$

The converse is just as important: if $p\ge 1$ and $d,e$ are positive divisors with $de=p^2+1$, then $(-p,p+d,p+e)$ satisfies the defining identity. So this is not merely a source of examples; it is the full classification exploited by the code.

Why only divisors up to $\sqrt{p^2+1}$ are needed

Once $p$ is fixed, the only freedom is the divisor pair $(d,e)$ of $p^2+1$. Swapping the two factors does not change the product, because

$$p(p+d)(p+e)=p(p+e)(p+d).$$

Hence it is enough to enumerate divisors $d\le \sqrt{p^2+1}$ and recover the partner $e=(p^2+1)/d$. This removes the trivial mirror duplication coming from $(d,e)$ and $(e,d)$.

Worked example: $p=3$

For $p=3$ we have

$$p^2+1=10.$$

The divisor pairs are $(1,10)$ and $(2,5)$. They produce the triples

$$(-3,4,13)\qquad\text{and}\qquad(-3,5,8),$$

so the corresponding Alexandrian integers are

$$3\cdot 4\cdot 13=156,\qquad 3\cdot 5\cdot 8=120.$$

Indeed,

$$-\frac{1}{3}+\frac{1}{4}+\frac{1}{13}=-\frac{1}{156},\qquad -\frac{1}{3}+\frac{1}{5}+\frac{1}{8}=-\frac{1}{120}.$$

This example explains two implementation choices. A single value of $p$ can generate several Alexandrian integers, and those values do not arrive in globally sorted order. The search therefore streams candidates into a structure that keeps only the smallest values seen so far.

A lower bound that makes termination safe

For every future parameter $p'$ and every admissible divisor pair, we always have $d\ge 1$ and $e\ge 1$. Therefore every unseen candidate obeys

$$A'=p'(p'+d)(p'+e)\ge p'(p'+1)^2.$$

After all candidates for a given $p$ have been processed, every later parameter $p+1,p+2,\dots$ satisfies

$$A'\ge (p+1)(p+2)^2.$$

If the container holding the current smallest $k$ candidates is already full and this lower bound is larger than its maximum element, then no future value can enter the smallest $k$. At that point the search may stop, and the current maximum is exactly the $k$-th Alexandrian integer.

How the Code Works

The C++, Python, and Java implementations iterate upward over $p=1,2,3,\dots$. For each $p$, they form $n=p^2+1$, enlarge an incremental prime table just enough to factor $n$, and then recursively generate all divisors $d\le \sqrt{n}$. Each such divisor determines the complementary divisor $e=n/d$ and therefore one candidate $A=p(p+d)(p+e)$.

Instead of storing every candidate ever seen, the implementations maintain a max-heap of fixed size $k$. While that heap is not yet full, every new unique candidate is inserted. Once it is full, only candidates smaller than the current heap maximum can matter; larger values are discarded immediately. A companion hash set tracks the values currently present in the heap so that equal Alexandrian integers are not counted twice.

The three versions differ mainly in numeric representation, not in mathematics. The C++ implementation uses 128-bit arithmetic for the products, the Python implementation relies on arbitrary-precision integers, and the Java implementation uses big integers. After each completed value of $p$, they apply the lower-bound test above. When that test proves that no later $p$ can improve the heap, the current heap maximum is returned.

Complexity Analysis

Let $P$ be the last parameter examined before the stopping criterion succeeds, let $\tau(n)$ be the divisor-counting function, and let $\pi(P)$ denote the number of stored primes up to $P$. For each $p\le P$, the algorithm factors $p^2+1$ by trial division through the current prime table and then enumerates the relevant divisors. A faithful high-level bound is

$$T(P)=O\!\left(\sum_{p=1}^{P}\bigl(\pi(p)+\tau(p^2+1)\log k\bigr)\right).$$

The $\pi(p)$ term reflects prime testing up to roughly $\sqrt{p^2+1}$, and the $\log k$ factor comes from possible heap updates when a candidate is admitted into the current best $k$.

Memory usage is

$$O\bigl(k+\pi(P)+\max_{p\le P}\tau(p^2+1)\bigr).$$

The dominant stored objects are the bounded heap, the hash set of active values, the growing prime table, and a temporary divisor buffer for the current factorization. The whole method is efficient because it replaces a three-variable search with a one-parameter scan, divisor generation, and aggressive early stopping.

Footnotes and References

Project Euler problem page: Problem 221 - Alexandrian Integers
Alexandrian integers: Wikipedia - Alexandrian integer
Divisors and divisor counting: Wikipedia - Divisor
Enumerating divisors from a prime factorization: cp-algorithms - Number of divisors / sum of divisors
Diophantine equations: Wikipedia - Diophantine equation

Problem 221 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <unordered_set>
#include <vector>
#include <cmath>

namespace {

using u64 = std::uint64_t;
using u32 = std::uint32_t;
using u128 = unsigned __int128;

struct Options {
    int target = 150000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }

    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--target=", options.target)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.target >= 1;
}

struct U128Hash {
    std::size_t operator()(const u128 v) const noexcept {
        const std::uint64_t lo = static_cast<std::uint64_t>(v);
        const std::uint64_t hi = static_cast<std::uint64_t>(v >> 64);
        return static_cast<std::size_t>(lo ^ (hi * 0x9e3779b97f4a7c15ULL));
    }
};

struct PrimeTable {
    std::vector<u32> primes{2U};
    u32 next_candidate = 3U;

    void ensure_up_to(const u64 limit) {
        while (primes.back() < limit) {
            bool is_prime = true;
            for (const u32 p : primes) {
                const u64 pu = static_cast<u64>(p);
                if (pu * pu > static_cast<u64>(next_candidate)) {
                    break;
                }
                if (next_candidate % p == 0U) {
                    is_prime = false;
                    break;
                }
            }
            if (is_prime) {
                primes.push_back(next_candidate);
            }
            next_candidate += 2U;
        }
    }
};

void build_divisors_leq(
    const std::vector<std::pair<u64, int>>& factors,
    const std::size_t idx,
    const u64 current,
    const u64 limit,
    std::vector<u64>& out
) {
    if (idx == factors.size()) {
        out.push_back(current);
        return;
    }

    const u64 prime = factors[idx].first;
    const int exponent = factors[idx].second;
    u64 power = 1ULL;

    for (int e = 0; e <= exponent; ++e) {
        if (current > limit / power) {
            break;
        }
        build_divisors_leq(factors, idx + 1U, current * power, limit, out);
        if (e == exponent || power > limit / prime) {
            break;
        }
        power *= prime;
    }
}

u64 nth_alexandrian(const int target) {
    std::vector<u128> heap;
    heap.reserve(static_cast<std::size_t>(target) + 16ULL);
    std::unordered_set<u128, U128Hash> active;
    active.reserve(static_cast<std::size_t>(target) * 2ULL + 64ULL);

    PrimeTable table;
    std::vector<std::pair<u64, int>> factors;
    factors.reserve(8U);
    std::vector<u64> divisors;
    divisors.reserve(64U);

    for (u64 p = 1ULL;; ++p) {
        const u64 n = p * p + 1ULL;
        const u64 root = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
        table.ensure_up_to(root + 1ULL);

        factors.clear();
        u64 rem = n;
        for (const u32 prime : table.primes) {
            const u64 q = static_cast<u64>(prime);
            if (q * q > rem) {
                break;
            }
            if (rem % q != 0ULL) {
                continue;
            }
            int count = 0;
            do {
                rem /= q;
                ++count;
            } while (rem % q == 0ULL);
            factors.push_back({q, count});
        }
        if (rem > 1ULL) {
            factors.push_back({rem, 1});
        }

        divisors.clear();
        build_divisors_leq(factors, 0U, 1ULL, root, divisors);

        for (const u64 d : divisors) {
            const u64 e = n / d;
            const u128 pu = static_cast<u128>(p);
            const u128 value = pu * (pu + static_cast<u128>(d)) * (pu + static_cast<u128>(e));

            if (heap.size() < static_cast<std::size_t>(target)) {
                if (active.insert(value).second) {
                    heap.push_back(value);
                    std::push_heap(heap.begin(), heap.end());
                }
                continue;
            }

            const u128 current_max = heap.front();
            if (value >= current_max) {
                continue;
            }
            if (active.find(value) != active.end()) {
                continue;
            }

            std::pop_heap(heap.begin(), heap.end());
            const u128 removed = heap.back();
            heap.pop_back();
            active.erase(removed);

            heap.push_back(value);
            std::push_heap(heap.begin(), heap.end());
            active.insert(value);
        }

        if (heap.size() == static_cast<std::size_t>(target)) {
            const u64 next_p = p + 1ULL;
            const u128 lower_bound = static_cast<u128>(next_p) *
                                     static_cast<u128>(next_p + 1ULL) *
                                     static_cast<u128>(next_p + 1ULL);
            if (lower_bound > heap.front()) {
                return static_cast<u64>(heap.front());
            }
        }
    }
}

bool run_checkpoints() {
    if (nth_alexandrian(6) != 630ULL) {
        std::cerr << "Checkpoint failed for 6th Alexandrian integer" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << nth_alexandrian(options.target) << '\n';
    return 0;
}

Python

import math
import heapq

def isqrt(n):
    if n < 0: return 0
    return int(math.sqrt(n))

def get_primes(limit):
    if limit < 2: return []
    is_prime = bytearray([1]) * (limit + 1)
    is_prime[0] = is_prime[1] = 0
    for p in range(2, isqrt(limit) + 1):
        if is_prime[p]:
            for i in range(p * p, limit + 1, p):
                is_prime[i] = 0
    return [p for p in range(2, limit + 1) if is_prime[p]]

def get_factors(n, primes):
    factors = []
    m = n
    for p in primes:
        if p * p > m: break
        if m % p == 0:
            count = 0
            while m % p == 0:
                m //= p
                count += 1
            factors.append((p, count))
    if m > 1:
        factors.append((m, 1))
    return factors

def get_divisors_leq(factors, idx, current, limit, out):
    if idx == len(factors):
        out.append(current)
        return
        
    p, exponent = factors[idx]
    power = 1
    for e in range(exponent + 1):
        if current > limit // power:
            break
        get_divisors_leq(factors, idx + 1, current * power, limit, out)
        if e == exponent or power > limit // p:
            break
        power *= p

def solve(target=150000):
    heap = [] # will store negative values to act as max-heap
    active = set()
    
    primes = [2, 3]
    
    def ensure_primes(limit):
        while primes[-1] < limit:
            cand = primes[-1] + 2
            while True:
                is_p = True
                for p in primes:
                    if p * p > cand: break
                    if cand % p == 0:
                        is_p = False
                        break
                if is_p:
                    primes.append(cand)
                    break
                cand += 2
                
    p = 1
    while True:
        n = p * p + 1
        root = isqrt(n)
        ensure_primes(root + 1)
        
        factors = get_factors(n, primes)
        
        divisors = []
        get_divisors_leq(factors, 0, 1, root, divisors)
        
        for d in divisors:
            e = n // d
            value = p * (p + d) * (p + e)
            
            if len(heap) < target:
                if value not in active:
                    heapq.heappush(heap, -value)
                    active.add(value)
                continue
                
            current_max = -heap[0]
            if value >= current_max:
                continue
            if value in active:
                continue
                
            removed = -heapq.heappop(heap)
            active.remove(removed)
            
            heapq.heappush(heap, -value)
            active.add(value)
            
        if len(heap) == target:
            next_p = p + 1
            lower_bound = next_p * (next_p + 1) * (next_p + 1)
            if lower_bound > -heap[0]:
                return str(-heap[0])
                
        p += 1

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;
import java.math.BigInteger;

public class Euler221 {
    static long isqrt(long n) {
        if (n < 0)
            return 0;
        long r = (long) Math.sqrt(n);
        while ((r + 1) * (r + 1) <= n)
            r++;
        while (r * r > n)
            r--;
        return r;
    }

    static class PrimeTable {
        List<Long> primes = new ArrayList<>();
        long nextCandidate = 3;

        PrimeTable() {
            primes.add(2L);
        }

        void ensureUpTo(long limit) {
            while (primes.get(primes.size() - 1) < limit) {
                boolean isPrime = true;
                for (long p : primes) {
                    if (p * p > nextCandidate)
                        break;
                    if (nextCandidate % p == 0) {
                        isPrime = false;
                        break;
                    }
                }
                if (isPrime)
                    primes.add(nextCandidate);
                nextCandidate += 2;
            }
        }
    }

    static class Factor {
        long p;
        int exp;

        Factor(long p, int exp) {
            this.p = p;
            this.exp = exp;
        }
    }

    static void buildDivisors(List<Factor> factors, int idx, long current, long limit, List<Long> out) {
        if (idx == factors.size()) {
            out.add(current);
            return;
        }

        long p = factors.get(idx).p;
        int exp = factors.get(idx).exp;
        long power = 1;

        for (int e = 0; e <= exp; ++e) {
            if (current > limit / power)
                break;
            buildDivisors(factors, idx + 1, current * power, limit, out);
            if (e == exp || power > limit / p)
                break;
            power *= p;
        }
    }

    public static String solve() {
        int target = 150000;
        PriorityQueue<BigInteger> heap = new PriorityQueue<>(Collections.reverseOrder());
        HashSet<BigInteger> active = new HashSet<>();
        PrimeTable table = new PrimeTable();

        long p = 1;
        while (true) {
            long n = p * p + 1;
            long root = isqrt(n);
            table.ensureUpTo(root + 1);

            List<Factor> factors = new ArrayList<>();
            long rem = n;
            for (long q : table.primes) {
                if (q * q > rem)
                    break;
                if (rem % q != 0)
                    continue;
                int count = 0;
                do {
                    rem /= q;
                    count++;
                } while (rem % q == 0);
                factors.add(new Factor(q, count));
            }
            if (rem > 1) {
                factors.add(new Factor(rem, 1));
            }

            List<Long> divisors = new ArrayList<>();
            buildDivisors(factors, 0, 1, root, divisors);

            for (long d : divisors) {
                long e = n / d;
                BigInteger bp = BigInteger.valueOf(p);
                BigInteger bd = BigInteger.valueOf(d);
                BigInteger be = BigInteger.valueOf(e);
                BigInteger value = bp.multiply(bp.add(bd)).multiply(bp.add(be));

                if (heap.size() < target) {
                    if (active.add(value)) {
                        heap.offer(value);
                    }
                    continue;
                }

                BigInteger currentMax = heap.peek();
                if (value.compareTo(currentMax) >= 0)
                    continue;
                if (active.contains(value))
                    continue;

                BigInteger removed = heap.poll();
                active.remove(removed);

                heap.offer(value);
                active.add(value);
            }

            if (heap.size() == target) {
                long nextP = p + 1;
                BigInteger bNext = BigInteger.valueOf(nextP);
                BigInteger bNextP1 = BigInteger.valueOf(nextP + 1);
                BigInteger lowerBound = bNext.multiply(bNextP1).multiply(bNextP1);

                if (lowerBound.compareTo(heap.peek()) > 0) {
                    return heap.peek().toString();
                }
            }

            p++;
        }
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}