Detailed mathematical approach

Problem Summary

Let \(P_n\) be the set of ordered pairs of lattice points \((a,b)\) present after \(n\) days, starting from the two initial states \(((1,0),(0,1))\) and \(((0,0),(0,0))\). The quantity asked for in the problem is \(g(n)=|P_n|\), and the implementations evaluate it at \(n=16\).

The naive process expands every state by extracting its first endpoint, second endpoint, and displacement vector, then rebuilding all states compatible with those three projections. That expansion is useful only for tiny \(n\). The real solution compresses one whole day into two scalar counts \(k_d\) and \(q_d\), then evaluates a closed formula in constant time.

Mathematical Approach

For a fixed day \(d\), write

\[ A_d=\{a:(a,b)\in P_d\},\qquad B_d=\{b:(a,b)\in P_d\},\qquad C_d=\{b-a:(a,b)\in P_d\}. \]

These are the left endpoints, right endpoints, and displacement vectors visible on day \(d\). The next day is built from exactly the three constructions used in the implementations.

The Three Derived Point Sets

If we denote the three day-\(d\) production rules by

\[ X_d=A_d\times B_d,\qquad Y_d=\{(b-c,b):b\in B_d,\ c\in C_d\},\qquad Z_d=\{(a,a+c):a\in A_d,\ c\in C_d\}, \]

then

\[ P_{d+1}=X_d\cup Y_d\cup Z_d. \]

Projecting this identity onto first coordinates, second coordinates, and differences gives three exact recurrences:

\[ A_{d+1}=B_d-C_d,\qquad B_{d+1}=A_d+C_d,\qquad C_{d+1}=B_d-A_d. \]

So the process is governed by three sets on the triangular lattice, not by an arbitrary cloud of pairs.

Why Only Two Numbers Matter on Each Day

The three projection sets are congruent by the symmetries of the update rule, so they always have the same size. Define

\[ k_d = |A_d| = |B_d| = |C_d|. \]

Each of \(X_d\), \(Y_d\), and \(Z_d\) therefore has size \(k_d^2\). The only remaining issue is overlap. A pair \((a,b)\) belongs to \(X_d\cap Y_d\) exactly when \(a\in A_d\), \(b\in B_d\), and \(b-a\in C_d\). But that same condition is also exactly what places \((a,b)\) in \(Z_d\). Hence all pairwise intersections and the triple intersection coincide. If we define

\[ q_d = \#\{(a,b): a\in A_d,\ b\in B_d,\ b-a\in C_d\}, \]

then inclusion-exclusion collapses to

\[ |P_{d+1}| = 3k_d^2 - 2q_d. \]

This is the origin of the formula used by all three implementations:

\[ g(0)=2,\qquad g(n)=3k_{n-1}^2-2q_{n-1}\quad (n\ge 1). \]

Hexagons on the Triangular Lattice

Once the recurrences for \(A_d,B_d,C_d\) are written out for a few small days, a rigid geometric pattern appears: each set is a lattice hexagon bounded by three strip constraints in the directions \(x\), \(y\), and \(x+y\). The three sets are the same hexagon seen under a \(120^\circ\) cyclic relabeling of those directions.

A convenient parameter is \(t\), which alternates with the day parity:

\[ t= \begin{cases} \dfrac{2^d-1}{3}, & d\ \text{even},\\[6pt] \dfrac{2^d-2}{3}, & d\ \text{odd}. \end{cases} \]

For example, on an even day the set \(A_d\) can be written as

\[ A_d=\{(x,y)\in\mathbb Z^2:\ -t\le x\le t+1,\ -t\le y\le t,\ -t\le x+y\le t+1\}, \]

and the other two sets are the same shape after cyclically shifting the three directions. On an odd day, the same description holds with the side-length triple changing from \((t,t,t+1)\) to \((t,t+1,t+1)\).

For a lattice hexagon with side-length triple \((u,v,w)\), the number of lattice points is

\[ uv+vw+wu+u+v+w+1. \]

Substituting \((u,v,w)=(t,t,t+1)\) or \((t,t+1,t+1)\) gives

\[ k_d= \begin{cases} 3t^2+5t+2, & d\ \text{even},\\[4pt] 3t^2+7t+4, & d\ \text{odd}. \end{cases} \]

Counting the Overlap Term \(q_d\)

The overlap term is a compatibility count between the three lattice hexagons. Fix \(c\in C_d\). Every valid pair \((a,b)\) with \(b-a=c\) is obtained by choosing \(a\in A_d\cap(B_d-c)\). Therefore

\[ q_d = \sum_{c\in C_d} |A_d\cap (B_d-c)|. \]

Because \(A_d\), \(B_d\), and \(C_d\) are aligned hexagons, each intersection \(A_d\cap(B_d-c)\) is again a smaller aligned hexagon or a boundary-degenerate version of one. Grouping the displacements \(c\) by their position in the hexagon turns the sum into quadratic row profiles, and summing those profiles yields the quartic closed forms

\[ q_d= \begin{cases} \dfrac{21t^4+70t^3+87t^2+46t+8}{4}, & d\ \text{even},\\[8pt] \dfrac{21t^4+98t^3+171t^2+134t+40}{4}, & d\ \text{odd}. \end{cases} \]

At that point the problem is finished: evaluate \(k_{n-1}\) and \(q_{n-1}\), then substitute them into \(g(n)=3k_{n-1}^2-2q_{n-1}\).

Worked Example: Day 1 to Day 2

After one day, the projection sets are

\[ A_1=\{(0,0),(0,1),(1,-1),(1,0)\}, \]

\[ B_1=\{(-1,1),(0,0),(0,1),(1,0)\}, \]

\[ C_1=\{(-1,0),(-1,1),(0,0),(0,1)\}. \]

So \(k_1=4\), and each production rule contributes \(4^2=16\) candidate pairs. The overlap term is

\[ q_1 = \sum_{c\in C_1} |A_1\cap(B_1-c)| = 2+3+3+2 = 10. \]

Therefore

\[ g(2)=|P_2| = 3\cdot 4^2 - 2\cdot 10 = 28, \]

which matches the small-value checks in the implementations. The same inclusion-exclusion mechanism is what drives the large-day formula.

How the Code Works

Brute-Force Validation for Tiny Days

The C++, Python, and Java implementations all use the same mathematics, but only the C++ version keeps a literal day-by-day set expansion for validation. It stores the current set \(P_d\), extracts \(A_d\), \(B_d\), and \(C_d\), applies the three production rules, and confirms that the closed formula reproduces the exact counts for the first few days.

Constant-Time Evaluation for the Real Input

The production path never materializes the lattice hexagons. It computes the day parameter \(d=n-1\), converts that to the parity-dependent value of \(t\), evaluates the polynomial formulas for \(k_d\) and \(q_d\), and finally returns \(3k_d^2-2q_d\). The C++ implementation uses 128-bit integers because the intermediate quartic expressions are already much larger than 64-bit signed limits; the Java version uses arbitrary-precision integers; Python gets exact big integers automatically.

All three implementations keep the easy base case \(g(0)=2\), and the C++ and Python implementations also verify small known values such as \(g(1)=8\) and \(g(2)=28\).

Complexity Analysis

The closed-form solver is \(O(1)\) time and \(O(1)\) memory: it evaluates a handful of powers, products, and additions on fixed-size formulas. The only reason the language implementations differ is numeric representation, not algorithmic structure.

The validation model is exponentially and then combinatorially explosive because it explicitly stores whole state sets \(P_d\). That version is intentionally confined to very small \(d\); it exists only to confirm the geometry and the inclusion-exclusion formula before the constant-time evaluator is used for the actual target \(n=16\).

Footnotes and References

1. Problem page: Project Euler 957: Point Genesis
2. Cartesian product: Wikipedia - Cartesian product
3. Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
4. Triangular lattice: Wikipedia - Triangular lattice
5. Lattice point: Wikipedia - Lattice point