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,$

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\).

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

If

$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

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