Problem 957: Point Genesis
View on Project EulerProject Euler Problem 957 Solution
EulerSolve provides an optimized solution for Project Euler Problem 957, Point Genesis, with C++, Python, Java, and a step-by-step mathematical explanation.
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....
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
- Problem page: Project Euler 957: Point Genesis
- Cartesian product: Wikipedia - Cartesian product
- Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
- Triangular lattice: Wikipedia - Triangular lattice
- Lattice point: Wikipedia - Lattice point
Problem 957 source code
C++
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <set>
#include <string>
#include <algorithm>
#include <functional>
namespace {
using i32 = std::int32_t;
using u128 = unsigned __int128;
u128 pow2_u128(int e) {
return (u128(1) << e);
}
std::string to_string_u128(u128 x) {
if (x == 0) {
return "0";
}
std::string s;
while (x > 0) {
const int d = static_cast<int>(x % 10);
s.push_back(static_cast<char>('0' + d));
x /= 10;
}
std::reverse(s.begin(), s.end());
return s;
}
std::pair<u128, u128> day_k_q(int day) {
if (day == 0) {
return {2, 2};
}
if ((day & 1) == 0) {
const u128 t = (pow2_u128(day) - 1) / 3;
const u128 t2 = t * t;
const u128 t3 = t2 * t;
const u128 t4 = t2 * t2;
const u128 k = 3 * t2 + 5 * t + 2;
const u128 q = (21 * t4 + 70 * t3 + 87 * t2 + 46 * t + 8) / 4;
return {k, q};
}
const u128 t = (pow2_u128(day) - 2) / 3;
const u128 t2 = t * t;
const u128 t3 = t2 * t;
const u128 t4 = t2 * t2;
const u128 k = 3 * t2 + 7 * t + 4;
const u128 q = (21 * t4 + 98 * t3 + 171 * t2 + 134 * t + 40) / 4;
return {k, q};
}
u128 g_formula(int n) {
if (n == 0) {
return 2;
}
const auto [k, q] = day_k_q(n - 1);
return 3 * k * k - 2 * q;
}
struct Vec2 {
i32 x{0};
i32 y{0};
bool operator<(const Vec2& other) const {
if (x != other.x) {
return x < other.x;
}
return y < other.y;
}
};
Vec2 add(Vec2 a, Vec2 b) {
return Vec2{static_cast<i32>(a.x + b.x), static_cast<i32>(a.y + b.y)};
}
Vec2 sub(Vec2 a, Vec2 b) {
return Vec2{static_cast<i32>(a.x - b.x), static_cast<i32>(a.y - b.y)};
}
struct PointState {
Vec2 a;
Vec2 b;
bool operator<(const PointState& other) const {
if (a.x != other.a.x) {
return a.x < other.a.x;
}
if (a.y != other.a.y) {
return a.y < other.a.y;
}
if (b.x != other.b.x) {
return b.x < other.b.x;
}
return b.y < other.b.y;
}
};
u128 g_bruteforce_small(int n) {
std::set<PointState> points;
points.insert(PointState{Vec2{1, 0}, Vec2{0, 1}});
points.insert(PointState{Vec2{0, 0}, Vec2{0, 0}});
for (int day = 0; day < n; ++day) {
std::set<Vec2> A;
std::set<Vec2> B;
std::set<Vec2> C;
for (const auto& p : points) {
A.insert(p.a);
B.insert(p.b);
C.insert(sub(p.b, p.a));
}
std::set<PointState> next;
for (const auto& a : A) {
for (const auto& b : B) {
next.insert(PointState{a, b});
}
}
for (const auto& b : B) {
for (const auto& c : C) {
next.insert(PointState{sub(b, c), b});
}
}
for (const auto& a : A) {
for (const auto& c : C) {
next.insert(PointState{a, add(a, c)});
}
}
points.swap(next);
}
return static_cast<u128>(points.size());
}
void run_validations() {
assert(g_formula(1) == 8);
assert(g_formula(2) == 28);
for (int n = 1; n <= 5; ++n) {
assert(g_formula(n) == g_bruteforce_small(n));
}
}
} // namespace
int main() {
run_validations();
std::cout << to_string_u128(g_formula(16)) << '\n';
return 0;
}
Python
def day_k_q(day):
if day == 0:
return 2, 2
if (day & 1) == 0:
t = ((1 << day) - 1) // 3
t2 = t * t
t3 = t2 * t
t4 = t2 * t2
k = 3 * t2 + 5 * t + 2
q = (21 * t4 + 70 * t3 + 87 * t2 + 46 * t + 8) // 4
return k, q
t = ((1 << day) - 2) // 3
t2 = t * t
t3 = t2 * t
t4 = t2 * t2
k = 3 * t2 + 7 * t + 4
q = (21 * t4 + 98 * t3 + 171 * t2 + 134 * t + 40) // 4
return k, q
def g_formula(n):
if n == 0:
return 2
k, q = day_k_q(n - 1)
return 3 * k * k - 2 * q
def solve():
return str(g_formula(16))
if __name__ == "__main__":
assert g_formula(1) == 8
assert g_formula(2) == 28
print(solve())
Java
import java.math.BigInteger;
public class Euler957 {
static BigInteger pow2(int e) {
return BigInteger.ONE.shiftLeft(e);
}
static BigInteger[] dayKQ(int day) {
if (day == 0)
return new BigInteger[] { BigInteger.TWO, BigInteger.TWO };
BigInteger TWO = BigInteger.TWO, THREE = BigInteger.valueOf(3), FIVE = BigInteger.valueOf(5);
BigInteger FOUR = BigInteger.valueOf(4), SEVEN = BigInteger.valueOf(7);
BigInteger EIGHT = BigInteger.valueOf(8);
BigInteger B21 = BigInteger.valueOf(21), B70 = BigInteger.valueOf(70);
BigInteger B87 = BigInteger.valueOf(87), B46 = BigInteger.valueOf(46);
BigInteger B98 = BigInteger.valueOf(98), B171 = BigInteger.valueOf(171);
BigInteger B134 = BigInteger.valueOf(134), B40 = BigInteger.valueOf(40);
if ((day & 1) == 0) {
BigInteger t = pow2(day).subtract(BigInteger.ONE).divide(THREE);
BigInteger t2 = t.multiply(t), t3 = t2.multiply(t), t4 = t2.multiply(t2);
BigInteger k = THREE.multiply(t2).add(FIVE.multiply(t)).add(TWO);
BigInteger q = B21.multiply(t4).add(B70.multiply(t3)).add(B87.multiply(t2)).add(B46.multiply(t)).add(EIGHT)
.divide(FOUR);
return new BigInteger[] { k, q };
} else {
BigInteger t = pow2(day).subtract(TWO).divide(THREE);
BigInteger t2 = t.multiply(t), t3 = t2.multiply(t), t4 = t2.multiply(t2);
BigInteger k = THREE.multiply(t2).add(SEVEN.multiply(t)).add(FOUR);
BigInteger q = B21.multiply(t4).add(B98.multiply(t3)).add(B171.multiply(t2)).add(B134.multiply(t)).add(B40)
.divide(FOUR);
return new BigInteger[] { k, q };
}
}
static BigInteger gFormula(int n) {
if (n == 0)
return BigInteger.TWO;
BigInteger[] kq = dayKQ(n - 1);
BigInteger k = kq[0], q = kq[1];
return BigInteger.valueOf(3).multiply(k).multiply(k).subtract(BigInteger.TWO.multiply(q));
}
public static void main(String[] args) {
System.out.println(gFormula(16));
}
}