Problem 604: Convex Path in Square
View on Project EulerProject Euler Problem 604 Solution
EulerSolve provides an optimized solution for Project Euler Problem 604, Convex Path in Square, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \(F(N)\) be the maximum number of lattice points inside an \(N \times N\) square that can lie on the graph of one strictly convex increasing function. If the graph passes through \(k+1\) lattice points, then those points determine \(k\) consecutive step vectors. The task is to evaluate \(F(10^{18})\), but the derivation below works for general \(N\). Mathematical Approach Write the selected lattice points in increasing \(x\)-order as \(P_0,P_1,\dots,P_k\). For each consecutive pair define the step vector $$P_i-P_{i-1}=(a_i,b_i).$$ Because the graph is increasing, every step satisfies \(a_i\gt 0\) and \(b_i\gt 0\). The entire problem becomes: how many such steps can we choose while preserving strict convexity and staying inside the square? Step 1: Convert the curve into step vectors Strict convexity means the slopes of consecutive chords are strictly increasing: $$\frac{b_1}{a_1} \lt \frac{b_2}{a_2} \lt \cdots \lt \frac{b_k}{a_k}.$$ The horizontal and vertical spans must both fit into the square, so $$\sum_{i=1}^{k} a_i \le N,\qquad \sum_{i=1}^{k} b_i \le N.$$ Adding the two inequalities gives the global \(L_1\) budget $$\sum_{i=1}^{k}(a_i+b_i)\le 2N.$$ Each extra segment increases the number of lattice points by exactly \(1\), so maximizing points is equivalent to maximizing the number of admissible step vectors....
Detailed mathematical approach
Problem Summary
Let \(F(N)\) be the maximum number of lattice points inside an \(N \times N\) square that can lie on the graph of one strictly convex increasing function. If the graph passes through \(k+1\) lattice points, then those points determine \(k\) consecutive step vectors. The task is to evaluate \(F(10^{18})\), but the derivation below works for general \(N\).
Mathematical Approach
Write the selected lattice points in increasing \(x\)-order as \(P_0,P_1,\dots,P_k\). For each consecutive pair define the step vector
$$P_i-P_{i-1}=(a_i,b_i).$$
Because the graph is increasing, every step satisfies \(a_i\gt 0\) and \(b_i\gt 0\). The entire problem becomes: how many such steps can we choose while preserving strict convexity and staying inside the square?
Step 1: Convert the curve into step vectors
Strict convexity means the slopes of consecutive chords are strictly increasing:
$$\frac{b_1}{a_1} \lt \frac{b_2}{a_2} \lt \cdots \lt \frac{b_k}{a_k}.$$
The horizontal and vertical spans must both fit into the square, so
$$\sum_{i=1}^{k} a_i \le N,\qquad \sum_{i=1}^{k} b_i \le N.$$
Adding the two inequalities gives the global \(L_1\) budget
$$\sum_{i=1}^{k}(a_i+b_i)\le 2N.$$
Each extra segment increases the number of lattice points by exactly \(1\), so maximizing points is equivalent to maximizing the number of admissible step vectors.
Step 2: Keep only primitive directions
If a step \((a,b)\) has \(d=\gcd(a,b)\gt 1\), then the reduced vector \((a/d,b/d)\) has the same slope but strictly smaller cost \((a+b)/d\). Since equal slopes cannot repeat on a strictly convex chain, a non-primitive step is never preferable to its primitive version. Therefore an optimal construction may be assumed to use only primitive vectors:
$$\gcd(a_i,b_i)=1.$$
Now fix
$$s=a+b.$$
Every primitive vector in that layer has the form \((a,s-a)\) with \(1\le a\le s-1\) and
$$\gcd(a,s)=1.$$
The number of such choices is exactly Euler's totient:
$$\#\{a:1\le a\le s-1,\ \gcd(a,s)=1\}=\varphi(s).$$
Step 3: Full layers are balanced
All primitive vectors with the same value of \(a+b=s\) form one layer. Every vector in that layer contributes one segment and costs \(s\) units of \(L_1\) budget.
The layer is symmetric because \((a,s-a)\) is primitive if and only if \((s-a,a)\) is primitive. Pairing those two vectors gives equal total contribution in the \(x\)- and \(y\)-directions, so
$$\sum_{\substack{1\le a\lt s\\ \gcd(a,s)=1}} a = \sum_{\substack{1\le a\lt s\\ \gcd(a,s)=1}} (s-a) = \frac{s\varphi(s)}{2}.$$
Therefore taking the whole layer uses total cost
$$s\varphi(s),$$
and that cost is split evenly between horizontal and vertical motion. This is the key reason the two-dimensional square constraint collapses to an \(L_1\)-budget problem for complete layers.
Step 4: Greedy by increasing \(a+b\)
Each available segment is worth one more lattice point, but a segment from layer \(s\) costs \(s\). Hence the cheapest segments always come from the smallest remaining \(s\), so the optimal count is obtained by taking layers in increasing order.
Define the prefix segment count and prefix cost by
$$\Phi(m)=\sum_{s=2}^{m}\varphi(s),\qquad C(m)=\sum_{s=2}^{m}s\varphi(s).$$
Let \(m\) be the largest integer such that
$$C(m)\le 2N.$$
Then all layers \(2,3,\dots,m\) are taken completely, contributing \(\Phi(m)\) segments. The next candidate layer is
$$s=m+1.$$
With remaining budget
$$R=2N-C(m),$$
the raw number of extra segments available from that next layer is
$$t=\left\lfloor\frac{R}{s}\right\rfloor,\qquad r=R\bmod s.$$
So the default segment count is \(\Phi(m)+t\), and the default point count is \(\Phi(m)+t+1\).
Step 5: Exact saturation is the only delicate case
When \(r\gt 0\), there is still slack after choosing \(t\) vectors from the final layer, and the implementations treat that situation as feasible. The only time a correction is needed is
$$r=0,\qquad t\equiv 1\pmod{2}.$$
In that case the lower complete layers are already perfectly balanced, and the last layer must balance itself exactly:
$$\sum a=\sum b=\frac{ts}{2}.$$
If \(t\) is even, we can use complementary pairs \((a,s-a)\) and \((s-a,a)\), so balancing is automatic.
If \(t=1\), balancing is impossible.
If \(s\equiv 0\pmod{4}\), balancing is also impossible. Every admissible \(a\) is odd, so the sum of an odd number of such terms is odd, while \(ts/2\) is even.
The remaining case is \(s\equiv 2\pmod{4}\), so \(s=2n\) with \(n\) odd. Any larger odd balanced subset can then be written as one balanced triple together with complementary pairs. So the question reduces to whether there exist three integers coprime to \(s\) whose sum is \(3s/2\). The implementation checks exactly that condition and subtracts one segment if it fails.
Worked Example: \(N=100\)
Here the total budget is
$$2N=200.$$
The prefix costs of the early layers are
$$\begin{aligned} C(2)&=2,\\ C(3)&=8,\\ C(4)&=16,\\ C(5)&=36,\\ C(6)&=48,\\ C(7)&=90,\\ C(8)&=122,\\ C(9)&=176. \end{aligned}$$
But \(C(10)=216\gt 200\), so the last complete layer is \(m=9\). The complete layers contribute
$$\Phi(9)=\varphi(2)+\varphi(3)+\cdots+\varphi(9)=1+2+2+4+2+6+4+6=27$$
segments. The remaining budget is
$$R=200-176=24.$$
The next layer is \(s=10\), hence
$$t=\left\lfloor\frac{24}{10}\right\rfloor=2,\qquad r=4.$$
No correction is needed because \(r\ne 0\). Therefore
$$F(100)=27+2+1=30,$$
matching the checkpoint used by the C++, Python, and Java implementations.
How the Code Works
The implementation first sieves Euler's totient function up to a cutoff large enough that the prefix cost \(C(m)\) exceeds \(2\times 10^{18}\). From those totients it builds two prefix arrays: one for \(\Phi(m)\) and one for \(C(m)\).
For a given \(N\), it binary-searches the largest \(m\) with \(C(m)\le 2N\). That immediately gives the number of complete layers and the remaining budget \(R\).
Next it computes \(t=\lfloor R/(m+1)\rfloor\) and \(r=R\bmod (m+1)\). The default answer is the number of complete-layer segments plus \(t\), followed by one extra point for the starting lattice point.
Finally, if the remainder is zero and the last-layer count is odd, the implementation runs the balancing test from Step 5. In the impossible subcases it decreases the segment count by \(1\), then returns the number of lattice points.
Complexity Analysis
If the sieve stops at \(L\), computing all totients up to \(L\) costs \(O(L\log\log L)\) time and \(O(L)\) memory. Building the prefix arrays is \(O(L)\). A query for \(F(N)\) then needs one binary search, a few arithmetic operations, and only a tiny extra check in the rare exact-saturation case.
Footnotes and References
- Problem page: https://projecteuler.net/problem=604
- Convex function: Wikipedia - Convex function
- Euler's totient function: Wikipedia - Euler's totient function
- Farey sequence: Wikipedia - Farey sequence
- Coprime integers: Wikipedia - Coprime integers
Problem 604 source code
C++
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>
// Project Euler 604: Maximum lattice points on a strictly convex increasing graph in an N x N square.
// After reducing to primitive step-vectors (a,b) with a,b>0 and gcd(a,b)=1, each step costs (a+b) in L1.
// Taking all steps with a+b=s gives exactly phi(s) choices and total L1 cost s*phi(s), so the L1-greedy
// optimum is found from prefix sums of s*phi(s), with a single parity/balancing exception when L1 is tight.
using u64 = std::uint64_t;
static bool has_balanced_triple_even2(u64 s) {
// s == 2*n with n odd. Need 3 totatives a,b,c of s summing to 3s/2.
const u64 n = s / 2;
if (n < 7) return false; // covers (6,10) style small failures quickly
const u64 LIM = 2000; // small even offsets are enough for n in our range
for (u64 d1 = 2; d1 <= LIM && d1 < n; d1 += 2) {
if (std::gcd(d1, n) != 1) continue;
for (u64 d2 = d1 + 2; d2 <= LIM && d1 + d2 < n; d2 += 2) {
if (std::gcd(d2, n) != 1) continue;
if (std::gcd(d1 + d2, n) != 1) continue;
return true;
}
}
return false;
}
static u64 F(u64 N, const std::vector<u64>& pref_phi, const std::vector<u64>& pref_cost) {
const u64 target = 2 * N; // total L1 budget
const auto it = std::upper_bound(pref_cost.begin(), pref_cost.end(), target);
const u64 m = (it == pref_cost.begin()) ? 0 : (u64)(it - pref_cost.begin() - 1);
const u64 base_cnt = (m >= 2) ? pref_phi[m] : 0; // sum_{s=2..m} phi(s)
const u64 base_cost = (m >= 2) ? pref_cost[m] : 0; // sum_{s=2..m} s*phi(s)
const u64 slack = target - base_cost;
const u64 s = m + 1; // next L1 layer
const u64 t = (s > 0) ? (slack / s) : 0;
const u64 r = (s > 0) ? (slack % s) : 0;
u64 segments = base_cnt + t;
// If r==0 and t is odd, L1 is saturated and we must have exact x=y=N; this requires a balanced odd subset
// in the last layer. For s divisible by 4 it's impossible by parity; for s==2*n (n odd) it's equivalent to
// existence of a balanced triple (then pad with complementary pairs).
if (r == 0 && (t & 1) && t > 0) {
bool ok = false;
if (t == 1) {
ok = false;
} else if ((s % 4) == 0) {
ok = false;
} else if ((s % 4) == 2) {
ok = has_balanced_triple_even2(s);
}
if (!ok) --segments;
}
return segments + 1; // lattice points = segments + start point
}
int main() {
const u64 N_max = 1000000000000000000ULL;
const u64 target_max = 2 * N_max;
// Sieve totients until prefix L1 cost covers 2*N_max.
u64 limit = 3000000;
std::vector<int> phi;
std::vector<u64> pref_phi;
std::vector<u64> pref_cost;
for (;;) {
phi.assign((size_t)limit + 1, 0);
for (u64 i = 0; i <= limit; ++i) phi[i] = (int)i;
for (u64 p = 2; p <= limit; ++p) {
if ((u64)phi[p] != p) continue;
for (u64 j = p; j <= limit; j += p) phi[j] -= phi[j] / (int)p;
}
pref_phi.assign((size_t)limit + 1, 0);
pref_cost.assign((size_t)limit + 1, 0);
for (u64 i = 2; i <= limit; ++i) {
pref_phi[i] = pref_phi[i - 1] + (u64)phi[i];
pref_cost[i] = pref_cost[i - 1] + i * (u64)phi[i];
}
if (pref_cost[limit] >= target_max) break;
limit = (u64)(limit * 1.4) + 1000;
}
// Statement validations.
if (F(1, pref_phi, pref_cost) != 2) {
std::cerr << "Validation failed: F(1)\n";
return 1;
}
if (F(3, pref_phi, pref_cost) != 3) {
std::cerr << "Validation failed: F(3)\n";
return 1;
}
if (F(9, pref_phi, pref_cost) != 6) {
std::cerr << "Validation failed: F(9)\n";
return 1;
}
if (F(11, pref_phi, pref_cost) != 7) {
std::cerr << "Validation failed: F(11)\n";
return 1;
}
if (F(100, pref_phi, pref_cost) != 30) {
std::cerr << "Validation failed: F(100)\n";
return 1;
}
if (F(50000, pref_phi, pref_cost) != 1898) {
std::cerr << "Validation failed: F(50000)\n";
return 1;
}
std::cout << F(N_max, pref_phi, pref_cost) << "\n";
return 0;
}
Python
import sys
import math
from bisect import bisect_right
def has_balanced_triple_even2(s):
n = s // 2
if n < 7: return False
LIM = 2000
for d1 in range(2, min(LIM + 1, n), 2):
if math.gcd(d1, n) != 1: continue
for d2 in range(d1 + 2, min(LIM + 1, n - d1), 2):
if math.gcd(d2, n) != 1: continue
if math.gcd(d1 + d2, n) != 1: continue
return True
return False
def F(N, pref_phi, pref_cost):
target = 2 * N
it = bisect_right(pref_cost, target)
m = 0 if it == 0 else it - 1
base_cnt = pref_phi[m] if m >= 2 else 0
base_cost = pref_cost[m] if m >= 2 else 0
slack = target - base_cost
s = m + 1
t = slack // s if s > 0 else 0
r = slack % s if s > 0 else 0
segments = base_cnt + t
if r == 0 and (t & 1) and t > 0:
ok = False
if t == 1:
ok = False
elif s % 4 == 0:
ok = False
elif s % 4 == 2:
ok = has_balanced_triple_even2(s)
if not ok:
segments -= 1
return segments + 1
def solve_N(N_max):
target_max = 2 * N_max
limit = 3000000
while True:
phi = list(range(limit + 1))
for p in range(2, limit + 1):
if phi[p] == p:
for j in range(p, limit + 1, p):
phi[j] -= phi[j] // p
pref_phi = [0] * (limit + 1)
pref_cost = [0] * (limit + 1)
for i in range(2, limit + 1):
pref_phi[i] = pref_phi[i - 1] + phi[i]
pref_cost[i] = pref_cost[i - 1] + i * phi[i]
if pref_cost[limit] >= target_max:
break
limit = int(limit * 1.4) + 1000
return str(F(N_max, pref_phi, pref_cost))
def solve():
return solve_N(1000000000000000000)
if __name__ == '__main__':
print(solve())
Java
public class Euler604 {
static long gcd(long a, long b) {
while (b != 0) {
long t = a % b;
a = b;
b = t;
}
return a;
}
static boolean hasBalancedTripleEven2(long s) {
long n = s / 2;
if (n < 7)
return false;
long LIM = 2000;
for (long d1 = 2; d1 <= LIM && d1 < n; d1 += 2) {
if (gcd(d1, n) != 1)
continue;
for (long d2 = d1 + 2; d2 <= LIM && d1 + d2 < n; d2 += 2) {
if (gcd(d2, n) != 1)
continue;
if (gcd(d1 + d2, n) != 1)
continue;
return true;
}
}
return false;
}
static int upperBound(long[] a, long key) {
int low = 0;
int high = a.length - 1;
int maxLessOrEqual = -1;
while (low <= high) {
int mid = (low + high) >>> 1;
if (a[mid] <= key) {
maxLessOrEqual = mid;
low = mid + 1;
} else {
high = mid - 1;
}
}
return maxLessOrEqual == -1 ? 0 : maxLessOrEqual + 1;
}
static long F(long N, long[] prefPhi, long[] prefCost) {
long target = 2 * N;
int it = upperBound(prefCost, target);
int m = (it == 0) ? 0 : it - 1;
long baseCnt = (m >= 2) ? prefPhi[m] : 0;
long baseCost = (m >= 2) ? prefCost[m] : 0;
long slack = target - baseCost;
long s = m + 1;
long t = (s > 0) ? (slack / s) : 0;
long r = (s > 0) ? (slack % s) : 0;
long segments = baseCnt + t;
if (r == 0 && (t % 2 != 0) && t > 0) {
boolean ok = false;
if (t == 1) {
ok = false;
} else if (s % 4 == 0) {
ok = false;
} else if (s % 4 == 2) {
ok = hasBalancedTripleEven2(s);
}
if (!ok)
segments--;
}
return segments + 1;
}
public static String solve() {
long N_max = 1000000000000000000L;
long target_max = 2 * N_max;
int limit = 3000000;
long[] prefPhi;
long[] prefCost;
while (true) {
int[] phi = new int[limit + 1];
for (int i = 0; i <= limit; i++)
phi[i] = i;
for (int p = 2; p <= limit; p++) {
if (phi[p] != p)
continue;
for (int j = p; j <= limit; j += p) {
phi[j] -= phi[j] / p;
}
}
prefPhi = new long[limit + 1];
prefCost = new long[limit + 1];
for (int i = 2; i <= limit; i++) {
prefPhi[i] = prefPhi[i - 1] + phi[i];
prefCost[i] = prefCost[i - 1] + (long) i * phi[i];
}
if (prefCost[limit] >= target_max)
break;
limit = (int) (limit * 1.4) + 1000;
}
return Long.toString(F(N_max, prefPhi, prefCost));
}
public static void main(String[] args) {
System.out.println(solve());
}
}