Problem 799: Pentagonal Puzzle

Project Euler Problem 799 Solution

EulerSolve provides an optimized solution for Project Euler Problem 799, Pentagonal Puzzle, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We seek the first pentagonal number $P_n=\frac{n(3n-1)}{2}$ that has more than 100 representations of the form $$P_n=P_a+P_b,\qquad a\ge b\ge 1.$$ The C++, Python, and Java implementations do not search over $(a,b)$ directly. Instead, they translate the question into a constrained sum-of-two-squares problem and then count the admissible representations through Gaussian-integer factorization. Mathematical Approach The starting point is the classical identity $$24P_t+1=(6t-1)^2.$$ That identity turns the pentagonal equation into an arithmetic problem about quadratic forms. Step 1: Convert the pentagonal equation into a sum of two squares If $P_a+P_b=P_n$, then $$24P_a+24P_b+2=24P_n+2,$$ so $$(6a-1)^2+(6b-1)^2=(6n-1)^2+1.$$ Define $$x=6a-1,\qquad y=6b-1,\qquad N_n=(6n-1)^2+1.$$ Then every pentagonal decomposition corresponds to a representation $$x^2+y^2=N_n$$ with $$x\equiv y\equiv 5\pmod 6,\qquad x\ge y\ge 5.$$ Conversely, every positive unordered pair $(x,y)$ satisfying these congruences maps back uniquely to $$a=\frac{x+1}{6},\qquad b=\frac{y+1}{6}.$$ So counting representations $P_n=P_a+P_b$ is exactly the same as counting unordered positive representations of $N_n$ as a sum of two squares with both coordinates congruent to $5$ modulo $6$....

Detailed mathematical approach

Problem Summary

We seek the first pentagonal number $P_n=\frac{n(3n-1)}{2}$ that has more than 100 representations of the form

$$P_n=P_a+P_b,\qquad a\ge b\ge 1.$$

The C++, Python, and Java implementations do not search over $(a,b)$ directly. Instead, they translate the question into a constrained sum-of-two-squares problem and then count the admissible representations through Gaussian-integer factorization.

Mathematical Approach

The starting point is the classical identity

$$24P_t+1=(6t-1)^2.$$

That identity turns the pentagonal equation into an arithmetic problem about quadratic forms.

Step 1: Convert the pentagonal equation into a sum of two squares

If $P_a+P_b=P_n$, then

$$24P_a+24P_b+2=24P_n+2,$$

$$(6a-1)^2+(6b-1)^2=(6n-1)^2+1.$$

Define

$$x=6a-1,\qquad y=6b-1,\qquad N_n=(6n-1)^2+1.$$

Then every pentagonal decomposition corresponds to a representation

$$x^2+y^2=N_n$$

with

$$x\equiv y\equiv 5\pmod 6,\qquad x\ge y\ge 5.$$

Conversely, every positive unordered pair $(x,y)$ satisfying these congruences maps back uniquely to

$$a=\frac{x+1}{6},\qquad b=\frac{y+1}{6}.$$

So counting representations $P_n=P_a+P_b$ is exactly the same as counting unordered positive representations of $N_n$ as a sum of two squares with both coordinates congruent to $5$ modulo $6$.

Step 2: Understand which prime factors can occur

Write

$$N_n=k^2+1,\qquad k=6n-1.$$

If an odd prime $q$ divides $N_n$, then

$$k^2\equiv -1\pmod q.$$

Therefore $-1$ is a quadratic residue modulo $q$, which can happen only when

$$q\equiv 1\pmod 4.$$

Also $k$ is odd, so $k^2\equiv 1\pmod 8$, hence

$$N_n=k^2+1\equiv 2\pmod 8.$$

This means that $2$ divides $N_n$ exactly once. Consequently every candidate that is fully factored has the shape

$$N_n=2\prod_{i=1}^{r} p_i^{e_i},\qquad p_i\equiv 1\pmod 4.$$

This is precisely the situation in which Gaussian integers give a natural parametrization of all representations by two squares.

Step 3: Generate the two-square representations in $\mathbb{Z}[i]$

For each odd prime $p_i\equiv 1\pmod 4$, choose integers $u_i,v_i>0$ such that

$$p_i=u_i^2+v_i^2=(u_i+iv_i)(u_i-iv_i).$$

Let $\pi_i=u_i+iv_i$. Since

$$2=(1+i)(1-i),$$

we can write

$$N_n=(1+i)(1-i)\prod_{i=1}^{r}\pi_i^{e_i}\overline{\pi_i}^{\,e_i}.$$

Now choose integers $s_i$ with $0\le s_i\le e_i$. The Gaussian integer

$$z=(1+i)\prod_{i=1}^{r}\pi_i^{s_i}\overline{\pi_i}^{\,e_i-s_i}$$

has norm

$$\mathcal{N}(z)=z\overline{z}=N_n.$$

If $z=x+iy$, then automatically

$$x^2+y^2=N_n.$$

The implementations enumerate exactly these exponent splits. Distinct primes contribute independently, so the total raw search space for one candidate is multiplicative in the exponents $e_i$.

Step 4: Filter to the pentagonal representations only

Not every representation of $N_n$ comes from pentagonal indices. The pair must satisfy the congruence condition

$$x\equiv y\equiv 5\pmod 6,$$

and both coordinates must be positive. Sign changes and swapping the coordinates do not create new decompositions of $P_n$, so the implementations convert each representation to the canonical ordered pair

$$\bigl(\max(|x|,|y|),\min(|x|,|y|)\bigr),$$

discard zero coordinates, and keep only distinct pairs that survive the modulo-$6$ test. Each surviving pair corresponds to exactly one choice of $(a,b)$ with $a\ge b$.

Step 5: Use a fast upper bound before exact enumeration

$$N_n=2\prod_{i=1}^{r} p_i^{e_i},$$

then the Gaussian construction above creates

$$2\prod_{i=1}^{r}(e_i+1)$$

signed and oriented products before deduplication, or equivalently $\prod (e_i+1)$ choices if the factor $2$ is included in the exponent list. Conjugate choices always collapse to the same unordered absolute pair, so the number of relevant unordered candidates is at most

$$\frac{1}{2}\prod_{j}(e_j+1),$$

where the product runs over all prime exponents, including the exponent $1$ of $2$. If this bound is at most $100$, then the current $n$ cannot solve the problem, so the exact Gaussian enumeration is skipped.

Worked Example: $n=49$

For $n=49$,

$$P_{49}=\frac{49(3\cdot 49-1)}{2}=3577,$$

and

$$N_{49}=(6\cdot 49-1)^2+1=293^2+1=85850=2\cdot 5^2\cdot 17\cdot 101.$$

The quick bound becomes

$$\frac{1}{2}(1+1)(2+1)(1+1)(1+1)=6,$$

so exact enumeration is necessary. After constructing all Gaussian products, removing sign and order duplicates, and imposing $x\equiv y\equiv 5\pmod 6$, the surviving pairs are

$$ (x,y)=(287,59)\quad\text{and}\quad (281,83). $$

They map back to

$$ (a,b)=\left(\frac{287+1}{6},\frac{59+1}{6}\right)=(48,10), $$

$$ (a,b)=\left(\frac{281+1}{6},\frac{83+1}{6}\right)=(47,14). $$

Hence

$$P_{49}=P_{48}+P_{10}=P_{47}+P_{14},$$

so $P_{49}$ has exactly two admissible pentagonal decompositions.

How the Code Works

The C++, Python, and Java implementations share the same core pipeline. They precompute a list of small primes and, for each prime $p\equiv 1\pmod 4$ in that table, one concrete decomposition $p=u^2+v^2$. Then they scan $n=1,2,3,\dots$, form $N_n=(6n-1)^2+1$, and trial-divide it by the precomputed primes. If a residual cofactor remains, the candidate is skipped; otherwise the full exponent list is known.

From that factorization, the implementation first evaluates the bound $\frac12\prod(e_j+1)$. Only if this can exceed $100$ does it build all Gaussian products coming from the exponent splits. Each product yields one pair $(x,y)$; the implementation then takes absolute values, orders the coordinates, removes zero coordinates, applies the congruence test $x\equiv y\equiv 5\pmod 6$, and stores the surviving pairs in a set. The first $n$ whose exact count is greater than $100$ produces the required answer $P_n$.

Complexity Analysis

For one candidate $n$, the cheap stage is trial division of $N_n$ by the precomputed prime table. The expensive stage, used only for candidates that survive the upper bound, enumerates $O\!\left(\prod(e_j+1)\right)$ Gaussian products and uses the same order of temporary storage before deduplication.

The total runtime depends on how far the search must advance before the first count above $100$ appears. In practice most candidates are rejected early: many are not fully factored by the fixed prime table, and many fully factored candidates fail the bound $\frac12\prod(e_j+1)>100$. That combination of cheap factor filtering and late exact counting is what makes the search feasible.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=799
Pentagonal numbers: Wikipedia - Pentagonal number
Fermat's theorem on sums of two squares: Wikipedia - Fermat's theorem on sums of two squares
Gaussian integers: Wikipedia - Gaussian integer
Quadratic residues: Wikipedia - Quadratic residue

Problem 799 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <map>
#include <set>
#include <utility>
#include <vector>

using u64 = std::uint64_t;
using i128 = __int128_t;

struct GaussInt {
    i128 x;
    i128 y;
};

static inline GaussInt multiply(const GaussInt& a, const GaussInt& b) {
    return {a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};
}

static inline GaussInt conjugate(const GaussInt& a) {
    return {a.x, -a.y};
}

static u64 isqrt_u64(const u64 n) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((r + 1ULL) * (r + 1ULL) <= n) ++r;
    while (r * r > n) --r;
    return r;
}

static std::vector<int> prime_list(const int limit) {
    std::vector<char> composite(static_cast<std::size_t>(limit + 1), 0);
    std::vector<int> primes;
    for (int i = 2; i <= limit; ++i) {
        if (composite[static_cast<std::size_t>(i)]) continue;
        primes.push_back(i);
        if (1LL * i * i <= limit) {
            for (int j = i * i; j <= limit; j += i) {
                composite[static_cast<std::size_t>(j)] = 1;
            }
        }
    }
    return primes;
}

static std::pair<int, int> sum_of_two_squares_prime(const int p) {
    for (int x = 1; 1LL * x * x <= p; ++x) {
        const int y2 = p - x * x;
        const u64 y = isqrt_u64(static_cast<u64>(y2));
        if (static_cast<int>(y * y) == y2) return {x, static_cast<int>(y)};
    }
    return {0, 0};
}

static u64 factor_with_primes(u64 value,
                              const std::vector<int>& primes,
                              const int max_prime,
                              std::vector<std::pair<int, int>>& factors) {
    factors.clear();
    for (const int p : primes) {
        if (p > max_prime) break;
        if (value % static_cast<u64>(p) != 0ULL) continue;
        int e = 0;
        do {
            value /= static_cast<u64>(p);
            ++e;
        } while (value % static_cast<u64>(p) == 0ULL);
        factors.push_back({p, e});
    }
    return value;
}

static int count_from_factorization(const std::vector<std::pair<int, int>>& factors,
                                    const std::map<int, GaussInt>& gaussian_primes) {
    std::vector<GaussInt> reps;
    reps.push_back({1, 0});

    for (const auto& [prime, exp] : factors) {
        const GaussInt base = gaussian_primes.at(prime);
        std::vector<GaussInt> powers(static_cast<std::size_t>(exp + 1));
        powers[0] = {1, 0};
        for (int i = 1; i <= exp; ++i) {
            powers[static_cast<std::size_t>(i)] = multiply(powers[static_cast<std::size_t>(i - 1)], base);
        }

        std::vector<GaussInt> terms(static_cast<std::size_t>(exp + 1));
        for (int i = 0; i <= exp; ++i) {
            terms[static_cast<std::size_t>(i)] =
                multiply(powers[static_cast<std::size_t>(i)],
                         conjugate(powers[static_cast<std::size_t>(exp - i)]));
        }

        std::vector<GaussInt> next;
        next.reserve(reps.size() * static_cast<std::size_t>(exp + 1));
        for (const auto& r : reps) {
            for (const auto& t : terms) {
                next.push_back(multiply(r, t));
            }
        }
        reps.swap(next);
    }

    std::set<std::pair<u64, u64>> uniq;
    for (const auto& rep : reps) {
        i128 x = rep.x < 0 ? -rep.x : rep.x;
        i128 y = rep.y < 0 ? -rep.y : rep.y;
        if (x == 0 || y == 0) continue;

        u64 xx = static_cast<u64>(x);
        u64 yy = static_cast<u64>(y);
        if (xx < yy) std::swap(xx, yy);
        if (xx % 6ULL == 5ULL && yy % 6ULL == 5ULL) {
            uniq.insert({xx, yy});
        }
    }
    return static_cast<int>(uniq.size());
}

static int brute_count(const int n) {
    auto pent = [](const int x) -> u64 { return static_cast<u64>(x) * (3ULL * static_cast<u64>(x) - 1ULL) / 2ULL; };
    const u64 target = pent(n);
    int cnt = 0;
    for (int a = 1; a <= n; ++a) {
        for (int b = 1; b <= a; ++b) {
            if (pent(a) + pent(b) == target) ++cnt;
        }
    }
    return cnt;
}

static int count_exact_small_n(const int n, std::map<int, GaussInt>& gaussian_primes) {
    const u64 k = 6ULL * static_cast<u64>(n) - 1ULL;
    u64 value = k * k + 1ULL;

    std::vector<std::pair<int, int>> factors;
    for (u64 p = 2; p * p <= value; ++p) {
        if (value % p != 0ULL) continue;
        int e = 0;
        do {
            value /= p;
            ++e;
        } while (value % p == 0ULL);

        const int ip = static_cast<int>(p);
        factors.push_back({ip, e});
        if (ip == 2) continue;
        if (ip % 4 == 3) return 0;
        if (!gaussian_primes.count(ip)) {
            const auto [x, y] = sum_of_two_squares_prime(ip);
            gaussian_primes[ip] = {x, y};
        }
    }
    if (value > 1ULL) {
        const int p = static_cast<int>(value);
        factors.push_back({p, 1});
        if (p != 2 && p % 4 == 3) return 0;
        if (p == 2) {
            gaussian_primes[p] = {1, 1};
        } else if (!gaussian_primes.count(p)) {
            const auto [x, y] = sum_of_two_squares_prime(p);
            gaussian_primes[p] = {x, y};
        }
    }
    return count_from_factorization(factors, gaussian_primes);
}

int main() {
    constexpr int threshold = 100;
    constexpr int prime_limit = 400;

    const std::vector<int> primes = prime_list(prime_limit);
    std::map<int, GaussInt> gaussian_primes;
    gaussian_primes[2] = {1, 1};

    for (const int p : primes) {
        if (p == 2 || p % 4 == 3) continue;
        const auto [x, y] = sum_of_two_squares_prime(p);
        gaussian_primes[p] = {x, y};
    }

    assert(count_exact_small_n(49, gaussian_primes) == 2);
    assert(count_exact_small_n(268, gaussian_primes) == 3);
    for (int n = 1; n <= 120; ++n) {
        assert(count_exact_small_n(n, gaussian_primes) == brute_count(n));
    }

    std::vector<std::pair<int, int>> factors;
    factors.reserve(16);

    for (u64 n = 1;; ++n) {
        const u64 k = 6ULL * n - 1ULL;
        const u64 value = k * k + 1ULL;

        const u64 rem = factor_with_primes(value, primes, prime_limit, factors);
        if (rem != 1ULL) continue;

        int bound = 1;
        for (const auto& [_, e] : factors) {
            bound *= (e + 1);
            if (bound > 2 * (threshold + 1)) break;
        }
        if (bound / 2 <= threshold) continue;

        const int ways = count_from_factorization(factors, gaussian_primes);
        if (ways > threshold) {
            const u64 answer = n * (3ULL * n - 1ULL) / 2ULL;
            std::cout << answer << '\n';
            return 0;
        }
    }
}

Python

import math

def solve():
    threshold = 100

    def isqrt(n):
        r = int(math.isqrt(n))
        while (r+1)*(r+1) <= n: r += 1
        while r*r > n: r -= 1
        return r

    def prime_list(limit):
        comp = bytearray(limit+1); ps = []
        for i in range(2, limit+1):
            if not comp[i]:
                ps.append(i)
                if i*i <= limit:
                    for j in range(i*i, limit+1, i): comp[j] = 1
        return ps

    def sum2sq(p):
        for x in range(1, int(p**0.5)+1):
            y2 = p - x*x; y = isqrt(y2)
            if y*y == y2: return (x, y)
        return (0, 0)

    def gm(a, b): return (a[0]*b[0]-a[1]*b[1], a[0]*b[1]+a[1]*b[0])
    def gconj(a): return (a[0], -a[1])

    primes = prime_list(400)
    gp = {2: (1,1)}
    for p in primes:
        if p == 2 or p % 4 == 3: continue
        gp[p] = sum2sq(p)

    def count_from_fac(factors):
        reps = [(1, 0)]
        for pr, exp in factors:
            base = gp[pr]
            pows = [(1, 0)]
            for i in range(exp): pows.append(gm(pows[-1], base))
            terms = [gm(pows[i], gconj(pows[exp-i])) for i in range(exp+1)]
            reps = [gm(r, t) for r in reps for t in terms]
        uniq = set()
        for x, y in reps:
            x, y = abs(x), abs(y)
            if x == 0 or y == 0: continue
            if x < y: x, y = y, x
            if x % 6 == 5 and y % 6 == 5: uniq.add((x, y))
        return len(uniq)

    for n in range(1, 10**9):
        k = 6*n - 1; value = k*k + 1
        factors = []
        v = value; ok = True
        for p in primes:
            if p > 400: break
            if v % p == 0:
                e = 0
                while v % p == 0: v //= p; e += 1
                factors.append((p, e))
                if p != 2 and p % 4 == 3: ok = False; break
                if p not in gp: gp[p] = sum2sq(p)
        if not ok: continue
        if v > 1: continue  # Only smooth numbers

        bound = 1
        for _, e in factors: bound *= (e+1)
        if bound // 2 <= threshold: continue

        ways = count_from_fac(factors)
        if ways > threshold:
            return str(n * (3*n - 1) // 2)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.HashSet;

public class Euler799 {

    static class GaussInt {
        BigInteger x, y;

        GaussInt(BigInteger x, BigInteger y) {
            this.x = x;
            this.y = y;
        }

        GaussInt(long x, long y) {
            this.x = BigInteger.valueOf(x);
            this.y = BigInteger.valueOf(y);
        }
    }

    static GaussInt multiply(GaussInt a, GaussInt b) {
        BigInteger axbx = a.x.multiply(b.x);
        BigInteger ayby = a.y.multiply(b.y);
        BigInteger axby = a.x.multiply(b.y);
        BigInteger aybx = a.y.multiply(b.x);
        return new GaussInt(axbx.subtract(ayby), axby.add(aybx));
    }

    static GaussInt conjugate(GaussInt a) {
        return new GaussInt(a.x, a.y.negate());
    }

    static ArrayList<Integer> primeList(int limit) {
        boolean[] composite = new boolean[limit + 1];
        ArrayList<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= limit; i++) {
            if (!composite[i]) {
                primes.add(i);
                if ((long) i * i <= limit) {
                    for (int j = i * i; j <= limit; j += i) {
                        composite[j] = true;
                    }
                }
            }
        }
        return primes;
    }

    static GaussInt sumOfTwoSquaresPrime(int p) {
        for (long x = 1; x * x <= p; x++) {
            long y2 = p - x * x;
            long y = (long) Math.sqrt(y2);
            if (y * y == y2) {
                return new GaussInt(x, y);
            }
        }
        return new GaussInt(0, 0);
    }

    static class Factor {
        int prime;
        int exp;

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

    static class Pair {
        long x, y;

        Pair(long x, long y) {
            this.x = x;
            this.y = y;
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (o == null || getClass() != o.getClass())
                return false;
            Pair pair = (Pair) o;
            return x == pair.x && y == pair.y;
        }

        @Override
        public int hashCode() {
            return java.util.Objects.hash(x, y);
        }
    }

    static long factorWithPrimes(long value, ArrayList<Integer> primes, int maxPrime, ArrayList<Factor> factors) {
        factors.clear();
        for (int p : primes) {
            if (p > maxPrime)
                break;
            if (value % p != 0)
                continue;
            int e = 0;
            do {
                value /= p;
                e++;
            } while (value % p == 0);
            factors.add(new Factor(p, e));
        }
        return value;
    }

    static int countFromFactorization(ArrayList<Factor> factors, HashMap<Integer, GaussInt> gaussianPrimes) {
        ArrayList<GaussInt> reps = new ArrayList<>();
        reps.add(new GaussInt(1, 0));

        for (Factor f : factors) {
            GaussInt base = gaussianPrimes.get(f.prime);
            GaussInt[] powers = new GaussInt[f.exp + 1];
            powers[0] = new GaussInt(1, 0);
            for (int i = 1; i <= f.exp; i++) {
                powers[i] = multiply(powers[i - 1], base);
            }

            GaussInt[] terms = new GaussInt[f.exp + 1];
            for (int i = 0; i <= f.exp; i++) {
                terms[i] = multiply(powers[i], conjugate(powers[f.exp - i]));
            }

            ArrayList<GaussInt> nextReps = new ArrayList<>();
            for (GaussInt r : reps) {
                for (GaussInt t : terms) {
                    nextReps.add(multiply(r, t));
                }
            }
            reps = nextReps;
        }

        HashSet<Pair> uniq = new HashSet<>();
        for (GaussInt rep : reps) {
            BigInteger bx = rep.x.abs();
            BigInteger by = rep.y.abs();

            if (bx.equals(BigInteger.ZERO) || by.equals(BigInteger.ZERO))
                continue;

            BigInteger bxx = bx.min(by);
            BigInteger byy = bx.max(by);

            BigInteger six = BigInteger.valueOf(6);
            BigInteger five = BigInteger.valueOf(5);

            if (bxx.remainder(six).equals(five) && byy.remainder(six).equals(five)) {
                uniq.add(new Pair(bxx.longValueExact(), byy.longValueExact()));
            }
        }
        return uniq.size();
    }

    static long solveN() {
        int threshold = 100;
        int primeLimit = 400;
        ArrayList<Integer> primes = primeList(primeLimit);

        HashMap<Integer, GaussInt> gaussianPrimes = new HashMap<>();
        gaussianPrimes.put(2, new GaussInt(1, 1));

        for (int p : primes) {
            if (p == 2 || p % 4 == 3)
                continue;
            gaussianPrimes.put(p, sumOfTwoSquaresPrime(p));
        }

        ArrayList<Factor> factors = new ArrayList<>();

        for (long n = 1;; n++) {
            long k = 6 * n - 1;
            long value = k * k + 1;

            long rem = factorWithPrimes(value, primes, primeLimit, factors);
            if (rem != 1)
                continue;

            int bound = 1;
            for (Factor f : factors) {
                bound *= (f.exp + 1);
                if (bound > 2 * (threshold + 1))
                    break;
            }
            if (bound / 2 <= threshold)
                continue;

            int ways = countFromFactorization(factors, gaussianPrimes);
            if (ways > threshold) {
                return n * (3 * n - 1) / 2;
            }
        }
    }

    public static String solve() {
        return Long.toString(solveN());
    }

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