Problem 566: Cake Icing Puzzle
View on Project EulerProject Euler Problem 566 Solution
EulerSolve provides an optimized solution for Project Euler Problem 566, Cake Icing Puzzle, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The puzzle is governed by three step sizes $$s_0=\frac{1}{a},\qquad s_1=\frac{1}{b},\qquad s_2=\frac{1}{\sqrt{c}},$$ applied cyclically to a point on the unit interval. For each triple \(9 \le a < b < c\), we seek the least positive number of steps \(T(a,b,c)\) for which the full three-phase system is synchronized again. The overall quantity is $$\sum_{9\le a < b < c \le n} T(a,b,c).$$ A brute-force simulation is hopeless because the synchronization time can be very large, so the solution converts the continuous process into a finite automaton and then solves the resulting congruence conditions exactly. Mathematical Approach The key observation is that each phase uses the same piecewise map with a different step size. For \(0<s<1\), define $$\Phi_s(x)=\begin{cases} 1-x, & 0\le x < s,\\ x-s, & s\le x < 1. \end{cases}$$ The dynamics applies \(\Phi_{s_0}\), then \(\Phi_{s_1}\), then \(\Phi_{s_2}\), and repeats this three-phase cycle forever. Each branch either reflects the interval orientation or preserves it, so every transition naturally carries one binary flip bit. Step 1: Turn the Continuous Process into a State Graph In each phase, the relevant breakpoints are the discontinuity \(s_\phi\), the endpoint \(0\), and all of their backward images under earlier phases....
Detailed mathematical approach
Problem Summary
The puzzle is governed by three step sizes
$$s_0=\frac{1}{a},\qquad s_1=\frac{1}{b},\qquad s_2=\frac{1}{\sqrt{c}},$$
applied cyclically to a point on the unit interval. For each triple \(9 \le a < b < c\), we seek the least positive number of steps \(T(a,b,c)\) for which the full three-phase system is synchronized again. The overall quantity is
$$\sum_{9\le a < b < c \le n} T(a,b,c).$$
A brute-force simulation is hopeless because the synchronization time can be very large, so the solution converts the continuous process into a finite automaton and then solves the resulting congruence conditions exactly.
Mathematical Approach
The key observation is that each phase uses the same piecewise map with a different step size. For \(0<s<1\), define
$$\Phi_s(x)=\begin{cases} 1-x, & 0\le x < s,\\ x-s, & s\le x < 1. \end{cases}$$
The dynamics applies \(\Phi_{s_0}\), then \(\Phi_{s_1}\), then \(\Phi_{s_2}\), and repeats this three-phase cycle forever. Each branch either reflects the interval orientation or preserves it, so every transition naturally carries one binary flip bit.
Step 1: Turn the Continuous Process into a State Graph
In each phase, the relevant breakpoints are the discontinuity \(s_\phi\), the endpoint \(0\), and all of their backward images under earlier phases. Once those breakpoints are known, every open interval between consecutive breakpoints behaves uniformly: every point in that interval goes to the same next interval and experiences the same orientation change.
Therefore one state is simply
$$\text{state}=(\text{phase},\ \text{interval}).$$
Because every state has exactly one successor, the entire system becomes a directed graph in which every component eventually lies on a cycle. The whole problem is then reduced to analyzing those cycles and determining when their orientation constraints can hold simultaneously.
Step 2: The Square Case Gives a Rational Grid
If \(c=k^2\), then \(s_2=1/k\) is rational. In that case all breakpoints lie on the common grid with denominator
$$D=\operatorname{lcm}(a,b,k).$$
Write
$$S_0=\frac{D}{a},\qquad S_1=\frac{D}{b},\qquad S_2=\frac{D}{k}.$$
Each phase is partitioned into the \(D\) cells
$$I_j=\left[\frac{j}{D},\frac{j+1}{D}\right),\qquad 0\le j<D.$$
For phase \(\phi\), the map is explicit:
$$I_j\longmapsto\begin{cases} I_{D-1-j}, & j<S_\phi,\\ I_{j-S_\phi}, & j\ge S_\phi. \end{cases}$$
The first branch is a reflection and contributes flip bit \(1\); the second branch is a translation and contributes flip bit \(0\). Thus the square case produces a finite automaton with exactly \(3D\) states.
Step 3: The Non-square Case Lives in a Quadratic Field
If \(c\) is not a square, then \(1/\sqrt{c}\) is irrational, so a simple rational grid no longer exists. The implementations instead represent every breakpoint exactly in the form
$$x=\frac{P}{ab}+\frac{Q\sqrt{c}}{c}\pmod{1},\qquad P,Q\in\mathbb{Z}.$$
This works because the three step sizes are
$$\frac{1}{a}=\frac{b}{ab},\qquad \frac{1}{b}=\frac{a}{ab},\qquad \frac{1}{\sqrt{c}}=\frac{\sqrt{c}}{c},$$
so every backward image remains in the same lattice generated by \(1/(ab)\) and \(\sqrt{c}/c\).
For each phase we begin with the seed set \(\{0,s_\phi\}\). If \(y\) is already a breakpoint in phase \(\phi\), then its preimages under the previous phase are
$$x=y+s_{\phi-1}\qquad\text{when } y<1-s_{\phi-1},$$
and
$$x=1-y\qquad\text{when } y>1-s_{\phi-1}.$$
Repeating this predecessor closure until no new points appear yields three finite breakpoint sets. They are then sorted exactly by comparing the sign of expressions of the form
$$A+B\sqrt{c},$$
which avoids floating-point ambiguity. Adjacent breakpoints again define intervals, so the non-square branch also becomes a finite deterministic automaton.
Step 4: Every Cycle Produces a Flip Word
Because phases advance \(0\to1\to2\to0\), every directed cycle has length \(3m\) for some \(m\). Along that cycle the automaton records a binary word
$$f_0,f_1,\dots,f_{3m-1},\qquad f_i\in\{0,1\},$$
where \(f_i=1\) means that step used the reflecting branch. Group the bits by phase:
$$\alpha_i=f_{3i},\qquad \beta_i=f_{3i+1},\qquad \gamma_i=f_{3i+2}\qquad(0\le i<m).$$
The orientation change over one full three-step round is
$$w_i=\alpha_i\oplus\beta_i\oplus\gamma_i.$$
If the total parity of the whole cycle is even, returning to the same place also returns to the same orientation after one geometric lap. If it is odd, one extra lap is needed to restore orientation. Therefore the natural modulus attached to that cycle is
$$P=\begin{cases} 3m, & \displaystyle\bigoplus_{i=0}^{m-1} w_i=0,\\ 6m, & \displaystyle\bigoplus_{i=0}^{m-1} w_i=1. \end{cases}$$
Step 5: Convert a Cycle into Admissible Residue Classes
Suppose the desired synchronization time has the form
$$N=3q+R,\qquad R\in\{0,1,2\}.$$
The partial phase offset \(R\) contributes its own correction word. For each \(i\), define
$$t_i^{(0)}=0,\qquad t_i^{(1)}=\alpha_i,\qquad t_i^{(2)}=\alpha_i\oplus\beta_i.$$
Then the corrected round word is
$$e_i^{(R)}=w_i\oplus t_i^{(R)}\oplus t_{i+1}^{(R)}.$$
A residue class is valid exactly when this corrected word is a cyclic rotation of \(w\) and the initial orientation parity is also consistent. The implementations find all such rotations with Knuth-Morris-Pratt on the doubled word \(w\,w\), then verify the start-parity condition with prefix xors. The result of one cycle is therefore a finite set
$$N\equiv r \pmod{P},\qquad r\in R_{\text{cycle}}.$$
Step 6: Merge All Cycle Constraints with CRT
Different connected components of the state graph impose independent congruence conditions on the same unknown \(N\). If one cycle allows residues \(R_1\) modulo \(P_1\) and another allows residues \(R_2\) modulo \(P_2\), then every compatible pair is merged by the generalized Chinese Remainder Theorem.
After processing all cycles, the surviving residue classes describe every synchronized time. The least positive one is exactly \(T(a,b,c)\).
Worked Example: The Square Triple \((10,14,16)\)
Here \(c=16=4^2\), so the rational-grid branch applies. We obtain
$$D=\operatorname{lcm}(10,14,4)=140,$$
hence
$$S_0=\frac{140}{10}=14,\qquad S_1=\frac{140}{14}=10,\qquad S_2=\frac{140}{4}=35.$$
So each phase has \(140\) cells and the automaton has \(420\) states in total. In phase \(0\), cells \(0\) through \(13\) reflect to cells \(139\) down to \(126\), while cells \(14\) through \(139\) translate to cells \(0\) through \(125\). The same rule is repeated with thresholds \(10\) and \(35\) in the next two phases.
Decomposing that automaton into directed cycles, extracting the admissible congruence classes from each cycle, and merging them with CRT yields the least positive synchronization time
$$T(10,14,16)=506,$$
which matches the checkpoint used by the implementations.
How the Code Works
The C++, Python, and Java implementations all follow the same mathematical plan. First they distinguish the square and non-square cases. In the square case they build the \(3D\)-state automaton directly from explicit cell formulas. In the non-square case they build exact breakpoint sets in quadratic-field coordinates, sort them with exact comparisons, and then form the corresponding interval automaton.
Next the implementation walks through every connected component, extracts its directed cycle, records the flip word, and converts that word into allowed residue classes by a combination of prefix-xor bookkeeping, rotation matching, and generalized CRT. The smallest positive residue surviving all cycle merges is the answer for that triple. Finally the outer loops sum these values for all \(9\le a<b<c\le n\), using arbitrary-precision integer arithmetic for the final total. The Python implementation delegates the heavy numerical work to the same compiled core, so all three languages share the same underlying algorithm.
Complexity Analysis
For one triple, let \(m\) be the number of intervals per phase after discretization. The resulting automaton has \(3m\) states. Building the square-case automaton is \(O(m)\). In the non-square case, predecessor closure and exact sorting dominate, giving roughly \(O(m \log m)\) work once the breakpoint set is known. Cycle extraction is linear in the number of states, and the word analysis for a cycle of length \(\ell\) is \(O(\ell)\) because both prefix-xor processing and KMP are linear.
The CRT stage depends on how many residue classes survive after each merge; in practice those sets stay modest, so the dominant cost is constructing and scanning the automaton. Memory usage is \(O(m)\) per triple. The overall computation for the final sum is the sum of these costs over all \(\binom{n-8}{3}\) triples, and the outer parameter loop is parallelized in the compiled implementations.
Footnotes and References
- Project Euler problem page: https://projecteuler.net/problem=566
- Chinese remainder theorem: Wikipedia - Chinese remainder theorem
- Knuth-Morris-Pratt algorithm: Wikipedia - Knuth-Morris-Pratt algorithm
- Quadratic field: Wikipedia - Quadratic field
- Symbolic dynamics: Wikipedia - Symbolic dynamics
Problem 566 source code
C++
#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <thread>
#include <unordered_set>
#include <vector>
#include <boost/multiprecision/cpp_int.hpp>
using boost::multiprecision::cpp_int;
static inline long long llgcd(long long a, long long b) {
while (b != 0) {
long long t = a % b;
a = b;
b = t;
}
return a < 0 ? -a : a;
}
static inline long long lllcm(long long a, long long b) {
return a / llgcd(a, b) * b;
}
struct Point {
long long P = 0;
long long Q = 0;
bool operator==(const Point &o) const { return P == o.P && Q == o.Q; }
};
struct PointHash {
size_t operator()(const Point &p) const noexcept {
uint64_t x = static_cast<uint64_t>(p.P) * 0x9e3779b97f4a7c15ULL;
uint64_t y = static_cast<uint64_t>(p.Q) + 0x9e3779b97f4a7c15ULL;
x ^= y + 0x9e3779b97f4a7c15ULL + (x << 6) + (x >> 2);
return static_cast<size_t>(x);
}
};
static inline int sign_AB_sqrtc(__int128 A, __int128 B, int c) {
if (B == 0) {
if (A == 0) return 0;
return (A > 0) ? 1 : -1;
}
if (A == 0) return (B > 0) ? 1 : -1;
if ((A > 0 && B > 0) || (A < 0 && B < 0)) {
return (A > 0) ? 1 : -1;
}
__int128 A2 = A * A;
__int128 B2c = B * B * static_cast<__int128>(c);
if (A < 0 && B > 0) {
if (B2c == A2) return 0;
return (B2c > A2) ? 1 : -1;
}
// A > 0, B < 0
if (B2c == A2) return 0;
return (B2c > A2) ? -1 : 1;
}
static inline int cmpPoints(const Point &a, const Point &b, long long ab, int c) {
long long dP = a.P - b.P;
long long dQ = a.Q - b.Q;
__int128 A = static_cast<__int128>(dP) * c;
__int128 B = static_cast<__int128>(dQ) * ab;
return sign_AB_sqrtc(A, B, c);
}
struct FieldContext {
int a = 0;
int b = 0;
int c = 0;
long long ab = 0;
long double sqrtc = 0.0L;
};
static inline Point canonical01(long long P, long long Q, const FieldContext &ctx) {
long double approx = static_cast<long double>(P) / static_cast<long double>(ctx.ab) +
static_cast<long double>(Q) * ctx.sqrtc / static_cast<long double>(ctx.c);
long long n = static_cast<long long>(std::floor(approx));
P -= n * ctx.ab;
auto less0 = [&](long long PP, long long QQ) {
__int128 A = static_cast<__int128>(PP) * ctx.c;
__int128 B = static_cast<__int128>(QQ) * ctx.ab;
return sign_AB_sqrtc(A, B, ctx.c) < 0;
};
auto ge1 = [&](long long PP, long long QQ) {
__int128 A = static_cast<__int128>(PP - ctx.ab) * ctx.c;
__int128 B = static_cast<__int128>(QQ) * ctx.ab;
return sign_AB_sqrtc(A, B, ctx.c) >= 0;
};
while (less0(P, Q)) P += ctx.ab;
while (ge1(P, Q)) P -= ctx.ab;
return {P, Q};
}
static inline Point addP(const Point &x, const Point &s, const FieldContext &ctx) {
return canonical01(x.P + s.P, x.Q + s.Q, ctx);
}
static inline Point subP(const Point &x, const Point &s, const FieldContext &ctx) {
return canonical01(x.P - s.P, x.Q - s.Q, ctx);
}
static inline Point oneMinus(const Point &x, const FieldContext &ctx) {
return canonical01(ctx.ab - x.P, -x.Q, ctx);
}
static std::vector<int> kmpPrefix(const std::vector<uint8_t> &pat) {
int n = static_cast<int>(pat.size());
std::vector<int> pi(n, 0);
for (int i = 1, j = 0; i < n; ++i) {
while (j > 0 && pat[i] != pat[j]) j = pi[j - 1];
if (pat[i] == pat[j]) ++j;
pi[i] = j;
}
return pi;
}
static std::vector<int> findRotations(const std::vector<uint8_t> &w,
const std::vector<uint8_t> &e) {
int n = static_cast<int>(w.size());
std::vector<uint8_t> text(2 * n);
for (int i = 0; i < n; ++i) {
text[i] = w[i];
text[i + n] = w[i];
}
auto pi = kmpPrefix(e);
std::vector<int> shifts;
for (int i = 0, j = 0; i < 2 * n - 1; ++i) {
while (j > 0 && text[i] != e[j]) j = pi[j - 1];
if (text[i] == e[j]) ++j;
if (j == n) {
int pos = i - n + 1;
if (pos < n) {
int shift = (n - pos) % n;
shifts.push_back(shift);
}
j = pi[j - 1];
}
}
return shifts;
}
static inline long long egcd(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long x1, y1;
long long g = egcd(b, a % b, x1, y1);
x = y1;
y = x1 - (a / b) * y1;
return g;
}
static inline bool crtPair(long long a1, long long m1, long long a2, long long m2,
long long &a, long long &m) {
long long g = llgcd(m1, m2);
long long diff = a2 - a1;
if (diff % g != 0) return false;
long long m1g = m1 / g;
long long m2g = m2 / g;
long long x, y;
egcd(m1g, m2g, x, y);
x %= m2g;
if (x < 0) x += m2g;
__int128 t = static_cast<__int128>(diff / g) * x % m2g;
if (t < 0) t += m2g;
__int128 l = static_cast<__int128>(m1) / g * m2;
__int128 res = static_cast<__int128>(a1) + static_cast<__int128>(m1) * t;
long long mod = static_cast<long long>(l);
long long ans = static_cast<long long>(res % l);
if (ans < 0) ans += mod;
a = ans;
m = mod;
return true;
}
static void mergeCRT(std::vector<long long> &sols, long long &M,
const std::vector<long long> &residues, long long P) {
std::vector<long long> next;
next.reserve(sols.size() * residues.size());
long long newM = lllcm(M, P);
for (long long s : sols) {
for (long long r : residues) {
long long a, m;
if (crtPair(s, M, r, P, a, m)) {
next.push_back(a);
}
}
}
std::sort(next.begin(), next.end());
next.erase(std::unique(next.begin(), next.end()), next.end());
sols.swap(next);
M = newM;
}
static std::pair<long long, std::vector<long long>>
cycleConstraint(const std::vector<uint8_t> &f) {
int L = static_cast<int>(f.size());
int n = L / 3;
uint8_t total_f = 0;
for (uint8_t b : f) total_f ^= b;
long long P = (total_f == 0) ? L : 2LL * L;
std::vector<uint8_t> a(n), b(n), c(n);
for (int i = 0; i < n; ++i) {
a[i] = f[3 * i];
b[i] = f[3 * i + 1];
c[i] = f[3 * i + 2];
}
std::vector<uint8_t> w(n);
uint8_t total_w = 0;
for (int i = 0; i < n; ++i) {
w[i] = a[i] ^ b[i] ^ c[i];
total_w ^= w[i];
}
std::vector<uint8_t> pref(n + 1, 0);
for (int i = 0; i < n; ++i) pref[i + 1] = pref[i] ^ w[i];
auto H0 = [&](int s) { return pref[s]; };
std::vector<long long> good;
for (int R = 0; R < 3; ++R) {
std::vector<uint8_t> t(n);
if (R == 0) {
std::fill(t.begin(), t.end(), 0);
} else if (R == 1) {
t = a;
} else {
for (int i = 0; i < n; ++i) t[i] = a[i] ^ b[i];
}
std::vector<uint8_t> dt(n);
for (int i = 0; i < n; ++i) dt[i] = t[i] ^ t[(i + 1) % n];
std::vector<uint8_t> e(n);
for (int i = 0; i < n; ++i) e[i] = w[i] ^ dt[i];
auto shifts = findRotations(w, e);
for (int sShift : shifts) {
for (int q = 0; q < 2; ++q) {
int mTrip = sShift + q * n;
uint8_t h0 = H0(sShift) ^ static_cast<uint8_t>(q * total_w);
uint8_t target = t[mTrip % n];
if (h0 != target) continue;
long long N = 3LL * mTrip + R;
long long Nmod = N % P;
good.push_back(Nmod);
}
}
}
std::sort(good.begin(), good.end());
good.erase(std::unique(good.begin(), good.end()), good.end());
return {P, good};
}
static long long F_square(int a, int b, int c) {
int k = static_cast<int>(std::llround(std::sqrt(static_cast<long double>(c))));
long long D = 1;
D = lllcm(D, a);
D = lllcm(D, b);
D = lllcm(D, k);
long long S[3] = {D / a, D / b, D / k};
int m = static_cast<int>(D);
int N = 3 * m;
std::vector<int> nxt(N);
std::vector<uint8_t> flip(N);
for (int ph = 0; ph < 3; ++ph) {
long long Sph = S[ph];
for (int j = 0; j < m; ++j) {
int j2;
uint8_t fbit;
if (j < Sph) {
j2 = m - 1 - j;
fbit = 1;
} else {
j2 = j - static_cast<int>(Sph);
fbit = 0;
}
int cur = ph * m + j;
int nxtph = (ph + 1) % 3;
nxt[cur] = nxtph * m + j2;
flip[cur] = fbit;
}
}
std::vector<char> vis(N, 0);
std::vector<std::pair<long long, std::vector<long long>>> constraints;
for (int s = 0; s < N; ++s) {
if (vis[s]) continue;
std::vector<int> cyc;
int cur = s;
while (!vis[cur]) {
vis[cur] = 1;
cyc.push_back(cur);
cur = nxt[cur];
}
int rot = 0;
for (int i = 0; i < static_cast<int>(cyc.size()); ++i) {
if (cyc[i] < m) { rot = i; break; }
}
std::rotate(cyc.begin(), cyc.begin() + rot, cyc.end());
std::vector<uint8_t> f(cyc.size());
for (int i = 0; i < static_cast<int>(cyc.size()); ++i) {
f[i] = flip[cyc[i]];
}
constraints.push_back(cycleConstraint(f));
}
std::sort(constraints.begin(), constraints.end(),
[](const auto &x, const auto &y) { return x.first < y.first; });
long long M = 1;
std::vector<long long> sols = {0};
for (auto &cc : constraints) {
mergeCRT(sols, M, cc.second, cc.first);
}
long long best = std::numeric_limits<long long>::max();
for (long long r : sols) {
if (r == 0) continue;
if (r < best) best = r;
}
if (best == std::numeric_limits<long long>::max()) best = M;
return best;
}
static int findIntervalIndex(const std::vector<Point> &pts, const Point &x,
const FieldContext &ctx) {
int lo = 0;
int hi = static_cast<int>(pts.size());
while (lo + 1 < hi) {
int mid = (lo + hi) / 2;
if (cmpPoints(pts[mid], x, ctx.ab, ctx.c) <= 0) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
static long long F_nonsquare(int a, int b, int c) {
FieldContext ctx{a, b, c, static_cast<long long>(a) * b, std::sqrt(static_cast<long double>(c))};
Point s0 = canonical01(ctx.ab / a, 0, ctx);
Point s1 = canonical01(ctx.ab / b, 0, ctx);
Point s2 = canonical01(0, 1, ctx);
Point s[3] = {s0, s1, s2};
Point thr[3] = {oneMinus(s0, ctx), oneMinus(s1, ctx), oneMinus(s2, ctx)};
std::array<std::unordered_set<Point, PointHash>, 3> sets;
std::array<std::vector<Point>, 3> stack;
for (int ph = 0; ph < 3; ++ph) {
sets[ph].reserve(1024);
sets[ph].insert(Point{0, 0});
sets[ph].insert(s[ph]);
stack[ph].push_back(Point{0, 0});
stack[ph].push_back(s[ph]);
}
while (true) {
bool progressed = false;
for (int ph = 0; ph < 3; ++ph) {
while (!stack[ph].empty()) {
progressed = true;
Point y = stack[ph].back();
stack[ph].pop_back();
int prev = (ph + 2) % 3;
int cmp = cmpPoints(y, thr[prev], ctx.ab, ctx.c);
if (cmp < 0) {
Point x = addP(y, s[prev], ctx);
if (sets[prev].insert(x).second) stack[prev].push_back(x);
}
if (cmp > 0) {
Point x = oneMinus(y, ctx);
if (sets[prev].insert(x).second) stack[prev].push_back(x);
}
}
}
if (!progressed) break;
}
std::vector<Point> pts[3];
int m = -1;
for (int ph = 0; ph < 3; ++ph) {
pts[ph].assign(sets[ph].begin(), sets[ph].end());
std::sort(pts[ph].begin(), pts[ph].end(),
[&](const Point &u, const Point &v) {
return cmpPoints(u, v, ctx.ab, ctx.c) < 0;
});
if (m < 0) m = static_cast<int>(pts[ph].size());
if (m != static_cast<int>(pts[ph].size())) {
std::cerr << "Partition size mismatch for (" << a << "," << b << "," << c << ")\n";
return -1;
}
}
int N = 3 * m;
std::vector<int> nxt(N);
std::vector<uint8_t> flip(N);
for (int ph = 0; ph < 3; ++ph) {
for (int i = 0; i < m; ++i) {
Point l = pts[ph][i];
Point r = (i + 1 < m) ? pts[ph][i + 1] : Point{ctx.ab, 0};
bool inside = cmpPoints(l, s[ph], ctx.ab, ctx.c) < 0;
Point yLeft = inside ? oneMinus(r, ctx) : subP(l, s[ph], ctx);
int nxtph = (ph + 1) % 3;
int j = findIntervalIndex(pts[nxtph], yLeft, ctx);
int cur = ph * m + i;
nxt[cur] = nxtph * m + j;
flip[cur] = inside ? 1 : 0;
}
}
std::vector<char> vis(N, 0);
std::vector<std::pair<long long, std::vector<long long>>> constraints;
for (int sIdx = 0; sIdx < N; ++sIdx) {
if (vis[sIdx]) continue;
std::vector<int> cyc;
int cur = sIdx;
while (!vis[cur]) {
vis[cur] = 1;
cyc.push_back(cur);
cur = nxt[cur];
}
int rot = 0;
for (int i = 0; i < static_cast<int>(cyc.size()); ++i) {
if (cyc[i] < m) { rot = i; break; }
}
std::rotate(cyc.begin(), cyc.begin() + rot, cyc.end());
std::vector<uint8_t> f(cyc.size());
for (int i = 0; i < static_cast<int>(cyc.size()); ++i) {
f[i] = flip[cyc[i]];
}
constraints.push_back(cycleConstraint(f));
}
std::sort(constraints.begin(), constraints.end(),
[](const auto &x, const auto &y) { return x.first < y.first; });
long long M = 1;
std::vector<long long> sols = {0};
for (auto &cc : constraints) {
mergeCRT(sols, M, cc.second, cc.first);
}
long long best = std::numeric_limits<long long>::max();
for (long long r : sols) {
if (r == 0) continue;
if (r < best) best = r;
}
if (best == std::numeric_limits<long long>::max()) best = M;
return best;
}
static inline bool isSquare(int c) {
int r = static_cast<int>(std::llround(std::sqrt(static_cast<long double>(c))));
return r * r == c;
}
static long long F(int a, int b, int c) {
if (isSquare(c)) return F_square(a, b, c);
return F_nonsquare(a, b, c);
}
static cpp_int computeG(int n, int threads) {
std::atomic<int> next_a(9);
std::vector<cpp_int> partial(threads);
std::vector<std::thread> workers;
workers.reserve(threads);
for (int t = 0; t < threads; ++t) {
workers.emplace_back([&, t]() {
cpp_int local = 0;
while (true) {
int a = next_a.fetch_add(1);
if (a > n - 2) break;
for (int b = a + 1; b <= n - 1; ++b) {
for (int c = b + 1; c <= n; ++c) {
local += F(a, b, c);
}
}
}
partial[t] = local;
});
}
for (auto &th : workers) th.join();
cpp_int sum = 0;
for (const auto &v : partial) sum += v;
return sum;
}
static bool validate_samples() {
struct SampleF { int a, b, c; long long expect; };
const SampleF f_samples[] = {
{9, 10, 11, 60},
{10, 14, 16, 506},
{15, 16, 17, 785232},
};
for (const auto &s : f_samples) {
long long got = F(s.a, s.b, s.c);
if (got != s.expect) {
std::cerr << "F(" << s.a << "," << s.b << "," << s.c << ") = " << got
<< ", expected " << s.expect << "\n";
return false;
}
}
struct SampleG { int n; long long expect; };
const SampleG g_samples[] = {
{11, 60},
{14, 58020},
{17, 1269260},
};
for (const auto &s : g_samples) {
cpp_int got = computeG(s.n, 1);
if (got != s.expect) {
std::cerr << "G(" << s.n << ") = " << got << ", expected " << s.expect << "\n";
return false;
}
}
return true;
}
int main() {
if (!validate_samples()) {
std::cerr << "Validation failed." << std::endl;
return 1;
}
int threads = static_cast<int>(std::thread::hardware_concurrency());
if (threads <= 0) threads = 1;
cpp_int ans = computeG(53, threads);
std::cout << ans << std::endl;
return 0;
}
Python
from __future__ import annotations
import re
import shutil
import subprocess
from pathlib import Path
ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")
def parse_output(stdout: str) -> str:
lines = [line.strip() for line in stdout.splitlines() if line.strip()]
if not lines:
return ""
answers = []
equals = []
for line in lines:
m1 = ANSWER_RE.search(line)
if m1:
answers.append(m1.group(1).strip())
m2 = EQUAL_RE.search(line)
if m2:
equals.append(m2.group(1).strip())
if answers:
return answers[-1]
if equals:
return equals[-1]
return lines[-1]
def should_skip_cpp_checkpoints(src: Path) -> bool:
try:
text = src.read_text(encoding="utf-8", errors="ignore")
except OSError:
return False
return "--skip-checkpoints" in text
def run_cpp(binary: Path, src: Path, root: Path) -> str:
cmd = [str(binary)]
if should_skip_cpp_checkpoints(src):
cmd.append("--skip-checkpoints")
try:
return subprocess.check_output(cmd, text=True, cwd=root)
except subprocess.CalledProcessError:
return subprocess.check_output(cmd, text=True, cwd=src.parent)
def solve() -> str:
problem_id = __file__.split("Euler")[-1].split(".")[0]
root = Path(__file__).resolve().parent.parent
src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"
if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
compiler = shutil.which("clang++") or shutil.which("g++")
if not compiler:
raise RuntimeError("No C++ compiler found (clang++/g++).")
subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])
output = run_cpp(binary=binary, src=src, root=root)
parsed = parse_output(output)
if not parsed:
raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
return parsed
if __name__ == "__main__":
print(solve())
Java
import java.math.BigInteger;
import java.util.*;
import java.util.concurrent.*;
public class Euler566 {
static long gcd(long a, long b) {
while (b != 0) {
long t = a % b;
a = b;
b = t;
}
return a < 0 ? -a : a;
}
static long lcm(long a, long b) {
if (a == 0 || b == 0)
return 0;
return (a / gcd(a, b)) * b;
}
static class Point {
long P, Q;
Point(long P, long Q) {
this.P = P;
this.Q = Q;
}
@Override
public boolean equals(Object o) {
if (this == o)
return true;
if (!(o instanceof Point))
return false;
Point point = (Point) o;
return P == point.P && Q == point.Q;
}
@Override
public int hashCode() {
long x = P * 0x9e3779b97f4a7c15L;
long y = Q + 0x9e3779b97f4a7c15L;
x ^= y + 0x9e3779b97f4a7c15L + (x << 6) + (x >>> 2);
return (int) (x ^ (x >>> 32));
}
}
static int signABsqrtc(long A, long B, int c) {
if (B == 0) {
if (A == 0)
return 0;
return (A > 0) ? 1 : -1;
}
if (A == 0)
return (B > 0) ? 1 : -1;
if ((A > 0 && B > 0) || (A < 0 && B < 0))
return (A > 0) ? 1 : -1;
long A2_lo = A * A;
long A2_hi = Math.multiplyHigh(A, A);
long B2_lo = B * B;
long B2_hi = Math.multiplyHigh(B, B);
long B2c_lo = B2_lo * c;
long B2c_hi = B2_hi * c + Math.multiplyHigh(B2_lo, c) + ((B2_lo < 0) ? c : 0);
int cmp;
if (B2c_hi != A2_hi) {
cmp = (B2c_hi > A2_hi) ? 1 : -1;
} else {
cmp = (B2c_lo + Long.MIN_VALUE > A2_lo + Long.MIN_VALUE) ? 1 : ((B2c_lo == A2_lo) ? 0 : -1);
}
if (A < 0 && B > 0) {
if (cmp == 0)
return 0;
return (cmp > 0) ? 1 : -1;
}
if (cmp == 0)
return 0;
return (cmp > 0) ? -1 : 1;
}
static int cmpPoints(Point a, Point b, long ab, int c) {
long dP = a.P - b.P;
long dQ = a.Q - b.Q;
long A = dP * c;
long B = dQ * ab;
return signABsqrtc(A, B, c);
}
static class FieldContext {
int a, b, c;
long ab;
double sqrtc;
FieldContext(int a, int b, int c) {
this.a = a;
this.b = b;
this.c = c;
this.ab = (long) a * b;
this.sqrtc = Math.sqrt(c);
}
}
static Point canonical01(long P, long Q, FieldContext ctx) {
double approx = ((double) P / ctx.ab) + ((double) Q * ctx.sqrtc / ctx.c);
long n = (long) Math.floor(approx);
P -= n * ctx.ab;
while (signABsqrtc(P * ctx.c, Q * ctx.ab, ctx.c) < 0) {
P += ctx.ab;
}
while (signABsqrtc((P - ctx.ab) * ctx.c, Q * ctx.ab, ctx.c) >= 0) {
P -= ctx.ab;
}
return new Point(P, Q);
}
static Point addP(Point x, Point s, FieldContext ctx) {
return canonical01(x.P + s.P, x.Q + s.Q, ctx);
}
static Point subP(Point x, Point s, FieldContext ctx) {
return canonical01(x.P - s.P, x.Q - s.Q, ctx);
}
static Point oneMinus(Point x, FieldContext ctx) {
return canonical01(ctx.ab - x.P, -x.Q, ctx);
}
static int[] kmpPrefix(byte[] pat) {
int n = pat.length;
int[] pi = new int[n];
for (int i = 1, j = 0; i < n; i++) {
while (j > 0 && pat[i] != pat[j])
j = pi[j - 1];
if (pat[i] == pat[j])
j++;
pi[i] = j;
}
return pi;
}
static List<Integer> findRotations(byte[] w, byte[] e) {
int n = w.length;
byte[] text = new byte[2 * n];
for (int i = 0; i < n; i++) {
text[i] = w[i];
text[i + n] = w[i];
}
int[] pi = kmpPrefix(e);
List<Integer> shifts = new ArrayList<>();
for (int i = 0, j = 0; i < 2 * n - 1; i++) {
while (j > 0 && text[i] != e[j])
j = pi[j - 1];
if (text[i] == e[j])
j++;
if (j == n) {
int pos = i - n + 1;
if (pos < n) {
shimsAdd(shifts, (n - pos) % n);
}
j = pi[j - 1];
}
}
return shifts;
}
static void shimsAdd(List<Integer> list, int val) {
list.add(val);
}
static long[] egcdExt(long a, long b) {
if (b == 0)
return new long[] { 1, 0, a };
long[] res = egcdExt(b, a % b);
long x1 = res[0];
long y1 = res[1];
long g = res[2];
long x = y1;
long y = x1 - (a / b) * y1;
return new long[] { x, y, g };
}
static long[] crtPair(long a1, long m1, long a2, long m2) {
long g = gcd(m1, m2);
long diff = a2 - a1;
if (diff % g != 0)
return null;
long m1g = m1 / g;
long m2g = m2 / g;
long[] res = egcdExt(m1g, m2g);
long x = res[0] % m2g;
if (x < 0)
x += m2g;
BigInteger bigX = BigInteger.valueOf(x);
BigInteger bigDiffG = BigInteger.valueOf(diff / g);
BigInteger bigM2g = BigInteger.valueOf(m2g);
BigInteger tBig = bigDiffG.multiply(bigX).remainder(bigM2g);
if (tBig.compareTo(BigInteger.ZERO) < 0)
tBig = tBig.add(bigM2g);
BigInteger m1G = BigInteger.valueOf(m1 / g);
BigInteger M = m1G.multiply(BigInteger.valueOf(m2));
BigInteger A = BigInteger.valueOf(a1).add(BigInteger.valueOf(m1).multiply(tBig)).remainder(M);
if (A.compareTo(BigInteger.ZERO) < 0)
A = A.add(M);
return new long[] { A.longValue(), M.longValue() };
}
static class CRTRes {
List<Long> sols;
long M;
CRTRes(List<Long> sols, long M) {
this.sols = sols;
this.M = M;
}
}
static CRTRes mergeCRT(List<Long> sols, long M, List<Long> residues, long P) {
List<Long> next = new ArrayList<>();
long newM = lcm(M, P);
for (long s : sols) {
for (long r : residues) {
long[] am = crtPair(s, M, r, P);
if (am != null) {
next.add(am[0]);
}
}
}
Set<Long> set = new HashSet<>(next);
next = new ArrayList<>(set);
Collections.sort(next);
return new CRTRes(next, newM);
}
static class CycleRes {
long P;
List<Long> good;
CycleRes(long P, List<Long> good) {
this.P = P;
this.good = good;
}
}
static CycleRes cycleConstraint(byte[] f) {
int L = f.length;
int n = L / 3;
byte totalF = 0;
for (byte b : f)
totalF ^= b;
long P = (totalF == 0) ? L : 2L * L;
byte[] a = new byte[n], b = new byte[n], c = new byte[n];
for (int i = 0; i < n; i++) {
a[i] = f[3 * i];
b[i] = f[3 * i + 1];
c[i] = f[3 * i + 2];
}
byte[] w = new byte[n];
byte totalW = 0;
for (int i = 0; i < n; i++) {
w[i] = (byte) (a[i] ^ b[i] ^ c[i]);
totalW ^= w[i];
}
byte[] pref = new byte[n + 1];
for (int i = 0; i < n; i++)
pref[i + 1] = (byte) (pref[i] ^ w[i]);
List<Long> good = new ArrayList<>();
for (int R = 0; R < 3; ++R) {
byte[] t = new byte[n];
if (R == 1) {
System.arraycopy(a, 0, t, 0, n);
} else if (R == 2) {
for (int i = 0; i < n; i++)
t[i] = (byte) (a[i] ^ b[i]);
}
byte[] dt = new byte[n];
for (int i = 0; i < n; i++)
dt[i] = (byte) (t[i] ^ t[(i + 1) % n]);
byte[] e = new byte[n];
for (int i = 0; i < n; i++)
e[i] = (byte) (w[i] ^ dt[i]);
List<Integer> shifts = findRotations(w, e);
for (int sShift : shifts) {
for (int q = 0; q < 2; ++q) {
int mTrip = sShift + q * n;
byte h0 = (byte) (pref[sShift] ^ (q * totalW % 2));
byte target = t[mTrip % n];
if (h0 != target)
continue;
long N_val = 3L * mTrip + R;
long Nmod = N_val % P;
good.add(Nmod);
}
}
}
Set<Long> set = new HashSet<>(good);
good = new ArrayList<>(set);
Collections.sort(good);
return new CycleRes(P, good);
}
static long F_square(int a, int b, int c) {
int k = (int) Math.round(Math.sqrt(c));
long D = 1;
D = lcm(D, a);
D = lcm(D, b);
D = lcm(D, k);
long[] S = { D / a, D / b, D / k };
int m = (int) D;
int N = 3 * m;
int[] nxt = new int[N];
byte[] flip = new byte[N];
for (int ph = 0; ph < 3; ++ph) {
long Sph = S[ph];
for (int j = 0; j < m; ++j) {
int j2;
byte fbit;
if (j < Sph) {
j2 = m - 1 - j;
fbit = 1;
} else {
j2 = j - (int) Sph;
fbit = 0;
}
int cur = ph * m + j;
int nxtph = (ph + 1) % 3;
nxt[cur] = nxtph * m + j2;
flip[cur] = fbit;
}
}
boolean[] vis = new boolean[N];
List<CycleRes> constraints = new ArrayList<>();
for (int s = 0; s < N; ++s) {
if (vis[s])
continue;
List<Integer> cyc = new ArrayList<>();
int cur = s;
while (!vis[cur]) {
vis[cur] = true;
cyc.add(cur);
cur = nxt[cur];
}
int rot = 0;
for (int i = 0; i < cyc.size(); i++) {
if (cyc.get(i) < m) {
rot = i;
break;
}
}
Collections.rotate(cyc, -rot);
byte[] fCyc = new byte[cyc.size()];
for (int i = 0; i < cyc.size(); i++)
fCyc[i] = flip[cyc.get(i)];
constraints.add(cycleConstraint(fCyc));
}
constraints.sort((x, y) -> Long.compare(x.P, y.P));
long M = 1;
List<Long> sols = new ArrayList<>();
sols.add(0L);
for (CycleRes cc : constraints) {
CRTRes r = mergeCRT(sols, M, cc.good, cc.P);
sols = r.sols;
M = r.M;
}
long best = Long.MAX_VALUE;
for (long r : sols) {
if (r == 0)
continue;
if (r < best)
best = r;
}
if (best == Long.MAX_VALUE)
best = M;
return best;
}
static int findIntervalIndex(List<Point> pts, Point x, FieldContext ctx) {
int lo = 0;
int hi = pts.size();
while (lo + 1 < hi) {
int mid = (lo + hi) / 2;
if (cmpPoints(pts.get(mid), x, ctx.ab, ctx.c) <= 0) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
static long F_nonsquare(int a, int b, int c) {
FieldContext ctx = new FieldContext(a, b, c);
Point s0 = canonical01(ctx.ab / a, 0, ctx);
Point s1 = canonical01(ctx.ab / b, 0, ctx);
Point s2 = canonical01(0, 1, ctx);
Point[] sArr = { s0, s1, s2 };
Point[] thr = { oneMinus(s0, ctx), oneMinus(s1, ctx), oneMinus(s2, ctx) };
@SuppressWarnings("unchecked")
Set<Point>[] sets = new HashSet[3];
@SuppressWarnings("unchecked")
List<Point>[] stack = new ArrayList[3];
for (int ph = 0; ph < 3; ++ph) {
sets[ph] = new HashSet<>();
stack[ph] = new ArrayList<>();
Point p0 = new Point(0, 0);
sets[ph].add(p0);
sets[ph].add(sArr[ph]);
stack[ph].add(p0);
stack[ph].add(sArr[ph]);
}
while (true) {
boolean progressed = false;
for (int ph = 0; ph < 3; ++ph) {
while (!stack[ph].isEmpty()) {
progressed = true;
Point y = stack[ph].remove(stack[ph].size() - 1);
int prev = (ph + 2) % 3;
int cmp = cmpPoints(y, thr[prev], ctx.ab, ctx.c);
if (cmp < 0) {
Point x = addP(y, sArr[prev], ctx);
if (sets[prev].add(x))
stack[prev].add(x);
}
if (cmp > 0) {
Point x = oneMinus(y, ctx);
if (sets[prev].add(x))
stack[prev].add(x);
}
}
}
if (!progressed)
break;
}
@SuppressWarnings("unchecked")
List<Point>[] pts = new ArrayList[3];
int m = -1;
for (int ph = 0; ph < 3; ++ph) {
pts[ph] = new ArrayList<>(sets[ph]);
pts[ph].sort((u, v) -> cmpPoints(u, v, ctx.ab, ctx.c));
if (m < 0)
m = pts[ph].size();
if (m != pts[ph].size()) {
return -1;
}
}
int N = 3 * m;
int[] nxt = new int[N];
byte[] flip = new byte[N];
for (int ph = 0; ph < 3; ++ph) {
for (int i = 0; i < m; ++i) {
Point l = pts[ph].get(i);
Point r = (i + 1 < m) ? pts[ph].get(i + 1) : new Point(ctx.ab, 0);
boolean inside = cmpPoints(l, sArr[ph], ctx.ab, ctx.c) < 0;
Point yLeft = inside ? oneMinus(r, ctx) : subP(l, sArr[ph], ctx);
int nxtph = (ph + 1) % 3;
int j = findIntervalIndex(pts[nxtph], yLeft, ctx);
int cur = ph * m + i;
nxt[cur] = nxtph * m + j;
flip[cur] = (byte) (inside ? 1 : 0);
}
}
boolean[] vis = new boolean[N];
List<CycleRes> constraints = new ArrayList<>();
for (int sIdx = 0; sIdx < N; ++sIdx) {
if (vis[sIdx])
continue;
List<Integer> cyc = new ArrayList<>();
int cur = sIdx;
while (!vis[cur]) {
vis[cur] = true;
cyc.add(cur);
cur = nxt[cur];
}
int rot = 0;
for (int i = 0; i < cyc.size(); i++) {
if (cyc.get(i) < m) {
rot = i;
break;
}
}
Collections.rotate(cyc, -rot);
byte[] fCyc = new byte[cyc.size()];
for (int i = 0; i < cyc.size(); i++)
fCyc[i] = flip[cyc.get(i)];
constraints.add(cycleConstraint(fCyc));
}
constraints.sort((x, y) -> Long.compare(x.P, y.P));
long M = 1;
List<Long> sols = new ArrayList<>();
sols.add(0L);
for (CycleRes cc : constraints) {
CRTRes r = mergeCRT(sols, M, cc.good, cc.P);
sols = r.sols;
M = r.M;
}
long best = Long.MAX_VALUE;
for (long r : sols) {
if (r == 0)
continue;
if (r < best)
best = r;
}
if (best == Long.MAX_VALUE)
best = M;
return best;
}
static boolean isSquare(int c) {
int r = (int) Math.round(Math.sqrt(c));
return r * r == c;
}
static long F(int a, int b, int c) {
if (isSquare(c))
return F_square(a, b, c);
return F_nonsquare(a, b, c);
}
static BigInteger computeG(int n) {
List<Callable<BigInteger>> tasks = new ArrayList<>();
for (int a = 9; a <= n - 2; a++) {
final int currA = a;
tasks.add(() -> {
BigInteger local = BigInteger.ZERO;
for (int b = currA + 1; b <= n - 1; b++) {
for (int c = b + 1; c <= n; c++) {
local = local.add(BigInteger.valueOf(F(currA, b, c)));
}
}
return local;
});
}
ExecutorService executor = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors());
BigInteger total = BigInteger.ZERO;
try {
for (Future<BigInteger> res : executor.invokeAll(tasks)) {
total = total.add(res.get());
}
} catch (Exception e) {
}
executor.shutdown();
return total;
}
public static String solve() {
return computeG(53).toString();
}
public static void main(String[] args) {
System.out.println(solve());
}
}