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

A positive integer $A$ is called an Alexandrian integer if there are integers $p$, $q$ and $r$ such that $A=pqr$ and

$$\frac1A=\frac1p+\frac1q+\frac1r.$$

For example, $630$ is one, with $p=5$, $q=-7$ and $r=-18$: indeed $5\cdot(-7)\cdot(-18)=630$ and $\frac15-\frac17-\frac1{18}=\frac1{630}$. Listed in increasing order, the first six Alexandrian integers are $6$, $42$, $120$, $156$, $420$ and $630$, so $630$ is the sixth. Each number appears once in this list, however many triples produce it, and we must find the $150000$th Alexandrian integer. ^[1]

The whole solution in one sentence: every Alexandrian integer equals $p(p+d)(p+e)$ for a positive integer $p$ and a factorisation $de=p^2+1$, so the program runs through $p=1,2,3,\dots$, factorises $p^2+1$, keeps the $150000$ smallest products found so far in a heap, and stops as soon as no larger $p$ can produce a product below the largest one kept.

Mathematical Approach

We turn the definition into a statement about the factorisations of $p^2+1$, check it on small cases, and then pick out the $k$-th smallest value, here $k=150000$, without generating all values.

1. Clearing the denominators

None of $p,q,r$ is $0$, because they are denominators. Multiplying $\frac1A=\frac1p+\frac1q+\frac1r$ by $A=pqr$ gives $1=qr+pr+pq$. So a positive integer $A$ is Alexandrian exactly when $A=pqr$ for integers $p,q,r$ with

$$pq+qr+rp=1.$$

For the example of the statement, $5\cdot(-7)+(-7)(-18)+(-18)\cdot5$${}=-35+126-90=1$.

2. Two of the three numbers are negative

If $p,q,r$ were all positive, each of the three products would be at least $1$ and their sum at least $3$. Since $A=pqr$ is positive, the number of negative factors is even, so exactly two of them are negative. Both conditions stay the same when $p,q,r$ are permuted, so we may rename them such that $p\gt0$ and $q,r\lt0$. Write $q=-u$ and $r=-v$ with positive integers $u$ and $v$; then $A=puv$, and the condition of Section 1 becomes $uv-pu-pv=1$.

3. A factorisation of $p^2+1$

Adding $p^2$ to both sides of $uv-pu-pv=1$ turns the left side into a product:

$$(u-p)(v-p)=p^2+1.$$

The two factors have the same sign, because their product is positive. They cannot both be negative: then $1\le p-u\le p-1$ and $1\le p-v\le p-1$, and their product would be at most $(p-1)^2\lt p^2+1$. Hence $d=u-p$ and $e=v-p$ are positive integers with $de=p^2+1$, and $A=p(p+d)(p+e)$.

Claim 1. A positive integer $A$ is an Alexandrian integer if and only if $A=p(p+d)(p+e)$ for positive integers $p,d,e$ with $de=p^2+1$.

Proof. Sections 1 to 3 prove one direction. Conversely, given such $p,d,e$, put $q=-(p+d)$ and $r=-(p+e)$. Then $pqr=p(p+d)(p+e)=A$ and

$$\begin{aligned}pq+qr+rp&=(p+d)(p+e)\\&\quad-p(2p+d+e)\\&=de-p^2=1,\end{aligned}$$

so $A$ is Alexandrian by Section 1.

The example of the statement is $p=5$ with $26=2\cdot13$, which gives $5\cdot7\cdot18=630$. The On-Line Encyclopedia of Integer Sequences describes the Alexandrian integers in the same way. ^[2]

4. The smallest cases

For $p=1,\dots,8$, the factorisations $de=p^2+1$ with $d\lt e$ and their products are:

$p=1$: $2=1\cdot2$ gives $1\cdot2\cdot3=6$;
$p=2$: $5=1\cdot5$ gives $2\cdot3\cdot7=42$;
$p=3$: $10=1\cdot10=2\cdot5$ gives $3\cdot4\cdot13=156$ and $3\cdot5\cdot8=120$;
$p=4$: $17=1\cdot17$ gives $4\cdot5\cdot21=420$;
$p=5$: $26=1\cdot26=2\cdot13$ gives $930$ and $630$;
$p=6$: $37=1\cdot37$ gives $6\cdot7\cdot43=1806$;
$p=7$: $50=1\cdot50=2\cdot25=5\cdot10$ gives $3192$, $2016$ and $7\cdot12\cdot17=1428$;
$p=8$: $65=1\cdot65=5\cdot13$ gives $5256$ and $8\cdot13\cdot21=2184$.

The products do not come out in increasing order: $p=7$ gives $1428$, less than $1806$ from $p=6$, and $p=8$ gives $2184$, less than $3192$ from $p=7$. Section 5 shows that every $p\ge9$ gives products of at least $3006$, so the Alexandrian integers below $3006$ are exactly $6$, $42$, $120$, $156$, $420$, $630$, $930$, $1428$, $1806$, $2016$ and $2184$, and the first six agree with the statement.

5. How large the product is for a given $p$

Since $de=p^2+1$, we have $(p+d)(p+e)$${}=p^2+p(d+e)+de$${}=2p^2+1+p(d+e)$, hence

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

For a fixed $p$ the product grows with $d+e$. Since $(d-1)(e-1)\ge0$, we have $d+e\le de+1=p^2+2$, with equality for $d=1$, which gives the largest product $A=p(p+1)(p^2+p+1)$. For the smallest, $(d+e)^2=(e-d)^2+4de$${}\ge4p^2+4\gt(2p)^2$, so $d+e\ge2p+1$ and

$$A\ge4p^3+p^2+p,$$

with equality when $d+e=2p+1$, as for $p=1,3,7$ with the products $6,120,1428$. For $p=9$ this bound is $3006$, and it increases with $p$. The program uses a simpler and weaker bound: $d\ge1$ and $e\ge1$ give $A\ge p(p+1)^2$, which also increases with $p$.

6. Only divisors up to $p$

Because $p^2\lt p^2+1\lt(p+1)^2$, the number $n=p^2+1$ is not a perfect square, so its divisors come in pairs $d\lt\sqrt n\lt e=n/d$, ^[3] and an integer below $\sqrt n$ is at most $p$. Exchanging $d$ and $e$ exchanges $q$ and $r$ and gives the same product. So each $p$ contributes one candidate for each divisor $d\le p$ of $n$, which makes $\tau(n)/2$ candidates, where $\tau(n)$ is the number of divisors of $n$. For $p=7$, the divisors of $50$ are $1,2,5,10,25,50$, and the three of them up to $7$ give the three products of Section 4.

7. Factorising $p^2+1$

Trial division ^[4] divides $n$ by the primes $2,3,5,\dots$ in turn, each as often as it divides, while the square of the prime does not exceed the remaining cofactor. The final cofactor is $1$ or a prime, since a composite cofactor $m$ has a prime factor at most $\sqrt m$, which would have been divided out already. Every cofactor is at most $p^2+1\lt(p+1)^2$, so the primes up to the first prime $\ge p+1$ suffice. The divisors of $n=\prod q_i^{a_i}$ are the products $\prod q_i^{b_i}$ with $0\le b_i\le a_i$. ^[3] When they are built prime by prime, a partial product above $p$ can be dropped at once, because further factors only enlarge it. For $p=7$ we get $50=2\cdot5^2$ and the divisors $1$, $5$ and $2$ up to $7$.

8. Keeping the $k$ smallest distinct values

A binary max-heap is a complete binary tree, filled level by level and stored in an array, in which every parent is at least as large as its children; the largest entry is at the root, and since a heap of $k$ entries has $\lfloor\log_2k\rfloor+1$ levels, inserting an entry or removing the root costs $O(\log k)$ steps. ^[5] The program keeps a max-heap $H$ of at most $k$ values and a hash set with the same values, and treats each candidate $A$ as follows. While $H$ has fewer than $k$ entries, $A$ is inserted unless it is already there. When $H$ is full, $A$ is rejected if it is at least the maximum $M$ of $H$ or already in $H$; otherwise $M$ is removed and $A$ is inserted.

Claim 2. After every step, $H$ consists of the $\min(k,s)$ smallest of the $s$ distinct values generated so far.

Proof. While $H$ is not full, nothing is removed, so $H$ holds every distinct value seen. Once it is full, its maximum $M$ never increases, because a removed maximum is replaced by a smaller value. Every earlier value that is not in $H$ is therefore at least $M$: it was rejected for being at least the maximum of its time, or it was removed as that maximum. A candidate above $M$ is not among the $k$ smallest, since $H$ holds $k$ smaller distinct values; a candidate equal to $M$ or found in $H$ is a repeat; and a candidate below $M$ that is not in $H$ is new, so replacing $M$ by it keeps exactly the $k$ smallest distinct values.

Claim 1 does not prove that different pairs $(p,d)$ give different products, so the set enforces distinctness without relying on such a proof. In fact each of the $1234824$ Alexandrian integers up to $10^{18}$ comes from a single pair $(p,d)$ with $d\le p$, and in the run for $k=150000$ the set test never rejects a candidate.

9. When to stop

Claim 3. If $H$ is full after all candidates of $p$ have been processed and $(p+1)(p+2)^2\gt M$, then $M$ is the $k$-th Alexandrian integer.

Proof. By Claim 1 and Section 6, every Alexandrian integer not generated yet comes from some $p'\ge p+1$, so by Section 5 it is at least $p'(p'+1)^2$${}\ge(p+1)(p+2)^2\gt M$. The $k$ smallest Alexandrian integers are therefore among the values generated so far, and by Claim 2 they are exactly the entries of $H$, the largest of which is $M$.

The test is eventually met: $H$ fills up, because the products $p(p+1)(p^2+p+1)$ for $d=1$ are all different; after that $M$ never increases, while $(p+1)(p+2)^2$ grows without bound.

10. Worked example: $k=6$

The candidates of $p\le4$ enter $H$ in the order $6,42,156,120,420$, and $930$ from $p=5$ fills it, with $M=930$. The second candidate of $p=5$, $630$, is smaller and not in $H$, so $930$ is removed and $630$ inserted; now $M=630$. The stopping test fails after $p=5$, because $6\cdot7^2=294$, and after $p=6$, whose candidate $1806$ is rejected, because $7\cdot8^2=448$. The three candidates $3192$, $1428$ and $2016$ of $p=7$ are rejected too, and now $8\cdot9^2=648\gt630$: the program stops and returns $630$, the sixth Alexandrian integer of the statement. With the bound $4p^3+p^2+p$ of Section 5, $4\cdot6^3+6^2+6=906\gt630$ would already end the search after $p=5$.

11. A link with formulas for $\pi$

Let $\alpha=\arctan\frac1{p+d}$ and $\beta=\arctan\frac1{p+e}$. By the addition formula for the tangent,

$$\begin{aligned}\tan(\alpha+\beta)&=\frac{\tan\alpha+\tan\beta}{1-\tan\alpha\tan\beta}\\&=\frac{2p+d+e}{(p+d)(p+e)-1},\end{aligned}$$

and the denominator $p^2+p(d+e)+de-1$ equals $p(2p+d+e)$, that is, $\tan(\alpha+\beta)=1/p$, exactly when $de=p^2+1$. Both angles lie between $0$ and $\pi/4$, so for positive integers $d,e$ the identity

$$\begin{aligned}\arctan\frac1p&=\arctan\frac1{p+d}\\&\quad+\arctan\frac1{p+e}\end{aligned}$$

holds exactly when $de=p^2+1$. For $p=1$ it is Euler's formula $\frac\pi4=\arctan\frac12+\arctan\frac13$, and $p=2$ gives $\arctan\frac12=\arctan\frac13+\arctan\frac17$; together they give $\frac\pi4=2\arctan\frac13$${}+\arctan\frac17$, another Machin-like formula for $\pi$. ^[6]

How the Code Works

The C++ program ^[7] starts in main: parse_arguments reads --target=, which sets $k$, $150000$ by default, and --skip-checkpoints; an unknown argument or $k\lt1$ gives exit code $1$, while a target above $2^{31}-1$ makes the int that accumulates its digits overflow, without any check. Unless skipped, run_checkpoints verifies that nth_alexandrian(6) returns $630$, with exit code $2$ on a failure. Then main prints nth_alexandrian(target), which carries out Sections 6 to 9 for $p=1,2,3,\dots$:

It sets $n=p^2+1$ and root to $\lfloor\sqrt n\rfloor=p$, computed in long double. PrimeTable::ensure_up_to(root + 1) extends the list of primes, testing the odd numbers in turn by trial division by the stored primes, until its last entry is at least $p+1$.
It factorises $n$ as in Section 7 into the pairs of prime and exponent in factors; a leftover cofactor above $1$ is a prime with exponent $1$.
build_divisors_leq chooses the exponent of each prime recursively and collects the divisors up to root. Its test current > limit / power detects that current * power exceeds the limit without forming the product, and cuts off that branch.
For each divisor $d$ it forms $e=n/d$ and the candidate $p(p+d)(p+e)$ as an unsigned __int128 and applies Section 8 with a std::vector kept in heap order by std::push_heap and std::pop_heap, with the maximum at the front, ^[8] and the std::unordered_set active, hashed by U128Hash.
After the last divisor it applies the test of Claim 3 and, if it succeeds, returns the maximum of the heap as a 64-bit integer.

The 128-bit products are needed: the $d=1$ candidate $p(p+1)(p^2+p+1)$ first exceeds $2^{64}$ at $p=2^{16}$, and the scan for $k=150000$ goes beyond that; with 64-bit products, $15$ wrapped-around values came out below the heap maximum in a test, and the printed value was wrong. The result itself fits into 64 bits: when the heap fills up, at $p=28037$, its entries are at most the $d=1$ product of that $p$, about $6.2\cdot10^{17}\lt2^{64}$, and the maximum never increases afterwards. ^[9]

Python and Java follow the same steps without the checkpoint. Python's solve(target=150000) takes the target as a parameter, which the script leaves at its default, while Java's solve() fixes it. Python's integers are unbounded; ^[9] its heapq min-heap holds the negated values, so the maximum is -heap[0], ^[8] its isqrt truncates math.sqrt, which gives $p$ exactly in the range scanned here, and its sieve get_primes is never called. Java does the bookkeeping in long and forms the candidates and the stopping bound as BigInteger values, kept in a PriorityQueue ordered by Collections.reverseOrder(), whose peek returns the maximum, and in a HashSet.

Complexity and Verification

Let $P$ be the last $p$ examined. Factorising $p^2+1$ tests at most $\pi(p)\approx p/\ln p$ primes, $\pi(x)$ being the number of primes up to $x$, and all of them when $p^2+1$ is prime; ^[4] repeated prime factors add only $O(\log p)$ divisions, since each at least halves the cofactor. So the factorisations cost $O(P^2/\log P)$ operations, and the prime table at most about $(P/2)\,\pi(\sqrt P)$ more. The divisors of $p^2+1$ take $O(\tau(p^2+1))$ steps, and each of the $D=\sum_{p\le P}\tau(p^2+1)/2$ candidates at most $O(\log k)$ heap steps and expected constant time in the set, which gives $O(P^2/\log P+D\log k)$ operations in all. For $k=150000$ we measured about $1.6\cdot10^8$ prime tests in the factorisations, $9\cdot10^5$ in the prime table and $D\approx7.5\cdot10^5$, so trial division dominates. The memory is $O(k)$ for the heap and the set, plus about $\pi(P)$ primes and the divisors of one $p$. On a desktop computer the C++ program took about $0.4$ seconds, Java about $1$ second and Python about $5$ seconds.

For $k=100$, $1000$ and $10000$ the scan ends after $p=105$, $916$ and $8566$, and there are $11684$, $121060$ and $1234824$ Alexandrian integers up to $10^{12}$, $10^{15}$ and $10^{18}$. These measurements suggest that $P$ grows about like $k$, but they prove no growth law, so the bounds above are stated in terms of $P$. The bound $4p^3+p^2+p$ of Section 5 would end the scan for $k=150000$ at about $63\%$ of the program's last $p$, with about $44\%$ of its prime tests. Since every odd prime factor of $p^2+1$ is $\equiv1\pmod 4$, ^[10] a sieve can also divide each such prime $q$ only out of the two residue classes of solutions of $x^2\equiv-1\pmod q$; the enumeration used for the checks below factorises $p^2+1$ in this way for all $p\lt10^6$ in well under a second.

The C++ program checks the sixth value $630$ of the statement on every run without --skip-checkpoints. For this article we also checked, outside the published programs:

a brute force over all integer triples with $|pqr|\le300000$ and every sign pattern with a positive product, using only the condition $pq+qr+rp=1$, finds exactly the $60$ Alexandrian integers up to $300000$ that Claim 1 gives, and no solution with three positive numbers;
a separate enumeration of all Alexandrian integers up to $10^{18}$, with a sieve factorisation of $p^2+1$ and a full sort instead of the heap and the stopping rule, gives the counts and the uniqueness quoted above; its 10th, 100th, 1000th and 10000th values $2016$, $1195116$, $772168278$ and $628833933090$ agree with the C++ program run with --target= and with the OEIS table, ^[2] and its $150000$th value with the output of the published programs;
the C++, Python and Java programs print the same value.

Footnotes and References

The references below supply the background facts used above; Claim 1, the bounds of Section 5 and Claims 2 and 3 are derived in this article.

Project Euler 221 — Alexandrian Integers. The official statement calls a positive integer A Alexandrian if A = pqr and 1/A = 1/p + 1/q + 1/r for some integers p, q, r, gives 630 (p = 5, q = −7, r = −18) as the sixth after 6, 42, 120, 156 and 420, and asks for the 150000th.
OEIS A147811 — Alexandrian integers. The entry defines them as the numbers n = pqr with 1/n = 1/p − 1/q − 1/r for integers p, q, r; a comment gives the form p(p + d)(p + (p² + 1)/d) with d a divisor of p² + 1, and a linked table lists the first 10000 terms.
Wikipedia — Divisor function. σ₀(n), also written τ(n), counts the divisors of n; it equals ∏ (aᵢ + 1) for n = ∏ pᵢ^aᵢ, and for a non-square n every divisor d is paired with n/d.
Wikipedia — Trial division. Trial factors need go no further than √n, because a composite number has a factor no larger than its square root.
Wikipedia — Binary heap. A binary heap is a complete binary tree, usually stored in an array; in a max-heap every key is at least as large as the keys of its children, and insertion and removal of the root take O(log n) time.
Wikipedia — Machin-like formula. Its two-term formulas include Euler's π/4 = arctan(1/2) + arctan(1/3) of 1737 and Hutton's or Vega's π/4 = 2 arctan(1/3) + arctan(1/7).
C++ — Euler221.cpp. The linked immutable C++ revision contains PrimeTable, build_divisors_leq, nth_alexandrian and run_checkpoints. It is the implementation record for this article.
cppreference — std::push_heap: with the default comparison the largest element is at the front. Python — heapq: heap[0] is the smallest item of a min-heap. Java — PriorityQueue: the head is the least element for the given ordering.
GCC manual — 128-bit Integers: unsigned __int128 is an unsigned 128-bit integer type. Python — numeric types: integers have unlimited precision. Java — BigInteger: immutable arbitrary-precision integers.
Wikipedia — Quadratic reciprocity. By the first supplement, for an odd prime p the congruence x² ≡ −1 (mod p) is solvable if and only if p ≡ 1 (mod 4).

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());
    }
}

Problem 221: Alexandrian Integers

Project Euler Problem 221 Solution

Detailed mathematical approach

Problem Summary

Mathematical Approach

1. Clearing the denominators

2. Two of the three numbers are negative

3. A factorisation of \(p^2+1\)

4. The smallest cases

5. How large the product is for a given \(p\)

6. Only divisors up to \(p\)

7. Factorising \(p^2+1\)

8. Keeping the \(k\) smallest distinct values