Problem 703: Circular Logic II
View on Project EulerProject Euler Problem 703 Solution
EulerSolve provides an optimized solution for Project Euler Problem 703, Circular Logic II, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For an \(n\)-bit state \(x\), write $$x=b_0+2b_1+4b_2+\cdots+2^{n-1}b_{n-1},\qquad b_i\in\{0,1\}.$$ The transition rule keeps shifting the bits to the right and inserts a new highest bit given by $$b_0\land(b_1\oplus b_2).$$ So the next state is $$f(x)=\left\lfloor\frac{x}{2}\right\rfloor+2^{n-1}\bigl(b_0\land(b_1\oplus b_2)\bigr).$$ The implementations reduce the counting task to the following graph problem on the \(N=2^n\) states: count all subsets of vertices such that no state is chosen together with its image under \(f\). In graph-theoretic language, this is an independent-set count on a special functional digraph, and the final result is taken modulo \(1001001011\). Mathematical Approach Let \(G_n\) be the directed graph with one vertex for each state \(x\in\{0,1,\dots,2^n-1\}\) and one directed edge \(x\to f(x)\). Because every vertex has exactly one outgoing edge, every connected component of \(G_n\) has a very rigid shape: one directed cycle, with rooted trees feeding into the cycle. Step 1: Turn the Bit Rule into a Functional Graph The mapping \(f\) is completely determined by the lowest three bits of \(x\). Every vertex therefore has outdegree \(1\), although its indegree may vary. This matters because graphs of outdegree \(1\) are functional graphs, and each component decomposes into exactly one cycle plus zero or more in-trees....
Detailed mathematical approach
Problem Summary
For an \(n\)-bit state \(x\), write
$$x=b_0+2b_1+4b_2+\cdots+2^{n-1}b_{n-1},\qquad b_i\in\{0,1\}.$$
The transition rule keeps shifting the bits to the right and inserts a new highest bit given by
$$b_0\land(b_1\oplus b_2).$$
So the next state is
$$f(x)=\left\lfloor\frac{x}{2}\right\rfloor+2^{n-1}\bigl(b_0\land(b_1\oplus b_2)\bigr).$$
The implementations reduce the counting task to the following graph problem on the \(N=2^n\) states: count all subsets of vertices such that no state is chosen together with its image under \(f\). In graph-theoretic language, this is an independent-set count on a special functional digraph, and the final result is taken modulo \(1001001011\).
Mathematical Approach
Let \(G_n\) be the directed graph with one vertex for each state \(x\in\{0,1,\dots,2^n-1\}\) and one directed edge \(x\to f(x)\). Because every vertex has exactly one outgoing edge, every connected component of \(G_n\) has a very rigid shape: one directed cycle, with rooted trees feeding into the cycle.
Step 1: Turn the Bit Rule into a Functional Graph
The mapping \(f\) is completely determined by the lowest three bits of \(x\). Every vertex therefore has outdegree \(1\), although its indegree may vary. This matters because graphs of outdegree \(1\) are functional graphs, and each component decomposes into exactly one cycle plus zero or more in-trees.
That decomposition is the structural reason the problem is tractable. Instead of counting admissible subsets on the full graph at once, we can solve each component separately and multiply the answers.
Step 2: Reformulate the Constraint as an Independent-Set Condition
If a subset \(S\) of states is valid, then no directed edge \(x\to f(x)\) may have both endpoints in \(S\). Since the restriction only concerns adjacent states, we can forget the arrow direction for counting purposes and work with the underlying undirected graph.
So the problem becomes: how many independent sets does each unicyclic component have? Once that is known, the total answer is
$$\prod_{\text{components }C}\operatorname{IS}(C)\pmod{1001001011},$$
where \(\operatorname{IS}(C)\) denotes the number of independent sets in component \(C\).
Step 3: Solve Every Tree Feeding into a Cycle
Take a non-cycle vertex \(u\), and orient its attached tree away from the cycle. Define
$$E(u)=\text{number of valid selections in the subtree of }u\text{ when }u\text{ is excluded},$$
$$I(u)=\text{number of valid selections in the subtree of }u\text{ when }u\text{ is included}.$$
If the children of \(u\) are \(v_1,\dots,v_k\), then standard tree DP gives
$$E(u)=\prod_{j=1}^{k}\bigl(E(v_j)+I(v_j)\bigr),$$
$$I(u)=\prod_{j=1}^{k}E(v_j).$$
The first formula allows each child to be either excluded or included. The second forbids including any child once \(u\) itself is included. A leaf satisfies
$$E(u)=1,\qquad I(u)=1.$$
Step 4: Collapse the Trees into Weights on the Cycle
Now consider a cycle \(c_1,c_2,\dots,c_m\). Each cycle vertex may have extra tree children that are not on the cycle. After solving those trees, we compress their effect into two local weights:
$$W_i^{(0)}=\prod_{v\in T_i}\bigl(E(v)+I(v)\bigr),$$
$$W_i^{(1)}=\prod_{v\in T_i}E(v),$$
where \(T_i\) is the set of non-cycle children attached directly to \(c_i\).
Interpretation:
\(W_i^{(0)}\) counts all valid choices in the incoming trees when \(c_i\) is not selected.
\(W_i^{(1)}\) counts all valid choices in the incoming trees when \(c_i\) is selected, so every attached child must stay out.
After this compression, the whole component is reduced to a weighted cycle.
Step 5: Run a Two-Case DP Around the Cycle
For a path version of the cycle, let
$$P_i^{(0)}=\text{weighted count up to }c_i\text{ with }c_i\text{ excluded},$$
$$P_i^{(1)}=\text{weighted count up to }c_i\text{ with }c_i\text{ included}.$$
Then for \(i\ge 2\),
$$P_i^{(0)}=\bigl(P_{i-1}^{(0)}+P_{i-1}^{(1)}\bigr)W_i^{(0)},$$
$$P_i^{(1)}=P_{i-1}^{(0)}W_i^{(1)}.$$
Because \(c_1\) and \(c_m\) are adjacent on the cycle, we split into two boundary cases.
Case A: \(c_1\) is excluded. Start with
$$P_1^{(0)}=W_1^{(0)},\qquad P_1^{(1)}=0.$$
This contributes
$$P_m^{(0)}+P_m^{(1)}.$$
Case B: \(c_1\) is included. Start with
$$P_1^{(0)}=0,\qquad P_1^{(1)}=W_1^{(1)}.$$
Now \(c_m\) must be excluded, so this contributes only
$$P_m^{(0)}.$$
If \(m=1\), the cycle is a self-loop. Then the lone cycle vertex cannot be selected at all, and the component contributes simply
$$W_1^{(0)}.$$
Worked Example: \(n=3\)
For \(n=3\) there are \(8\) states. Writing states in binary, the transition graph splits into two components:
000 -> 000, with the chain 100 -> 010 -> 001 -> 000.
011 -> 101 -> 110 -> 011, with the extra leaf 111 -> 011.
In the first component, 000 is a self-loop, so it cannot be chosen. What remains is the path 100-010-001, whose independent sets are
$$\varnothing,\ \{100\},\ \{010\},\ \{001\},\ \{100,001\},$$
so that component contributes \(5\).
In the second component, the leaf attached to 011 gives weights
$$W_1^{(0)}=2,\qquad W_1^{(1)}=1,$$
while the other two cycle vertices have weights
$$W_2^{(0)}=W_2^{(1)}=W_3^{(0)}=W_3^{(1)}=1.$$
Case A, first cycle vertex excluded:
$$P_1=(2,0),\qquad P_2=(2,2),\qquad P_3=(4,2),$$
so the contribution is \(4+2=6\).
Case B, first cycle vertex included:
$$P_1=(0,1),\qquad P_2=(1,0),\qquad P_3=(1,1),$$
and only the excluded-last state is allowed, so the contribution is \(1\).
Thus the 3-cycle component contributes \(6+1=7\), and the total is
$$5\cdot 7=35,$$
which matches the small checkpoint used by the implementation.
How the Code Works
The C++, Python, and Java implementations first build the transition table for all \(2^n\) states and, at the same time, the reverse adjacency lists that list every predecessor of each state. They then find every directed cycle by a standard three-color walk: unvisited, currently on the path, and fully processed. Whenever a walk re-enters a currently active vertex, the cycle segment is recorded.
After the cycles are known, the implementation evaluates every incoming tree exactly once with a post-order dynamic program. Those tree values are folded into the two weights for each cycle vertex, and then a second DP is run around the cycle using the two boundary cases described above. Each component count is multiplied into the global answer modulo \(1001001011\).
Complexity Analysis
Let \(N=2^n\). Building the transition table and reverse graph takes \(O(N)\) time and \(O(N)\) memory. Cycle detection also visits each vertex a constant number of times, so it is \(O(N)\). The tree DP and cycle DP together process each edge and each vertex only once more, so the total running time remains
$$O(N)=O(2^n),$$
with memory usage
$$O(N)=O(2^n).$$
For the actual input \(n=20\), this means linear work over \(1{,}048{,}576\) states, which is exactly why the graph decomposition is effective.
Footnotes and References
- Problem page: https://projecteuler.net/problem=703
- Independent set: Wikipedia - Independent set
- Dynamic programming: Wikipedia - Dynamic programming
- Cycle finding in directed graphs: CP-Algorithms - Checking a graph for acyclicity and finding a cycle
- Shift registers and feedback rules: Wikipedia - Shift register
Problem 703 source code
C++
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <utility>
#include <vector>
namespace {
using u64 = std::uint64_t;
constexpr u64 kMod = 1'001'001'011ULL;
u64 add_mod(const u64 a, const u64 b) {
const u64 s = a + b;
if (s >= kMod) {
return s - kMod;
}
return s;
}
u64 mul_mod(const u64 a, const u64 b) {
return static_cast<u64>((static_cast<unsigned __int128>(a) * static_cast<unsigned __int128>(b)) %
static_cast<unsigned __int128>(kMod));
}
u64 solve(const int n) {
const int N = 1 << n;
std::vector<int> succ(static_cast<std::size_t>(N), 0);
for (int x = 0; x < N; ++x) {
const int b1 = x & 1;
const int b2 = (x >> 1) & 1;
const int b3 = (x >> 2) & 1;
const int last = b1 & (b2 ^ b3);
succ[static_cast<std::size_t>(x)] = (x >> 1) | (last << (n - 1));
}
std::vector<int> head(static_cast<std::size_t>(N), -1);
std::vector<int> to(static_cast<std::size_t>(N), 0);
std::vector<int> next_edge(static_cast<std::size_t>(N), -1);
int edge_ptr = 0;
for (int v = 0; v < N; ++v) {
const int u = succ[static_cast<std::size_t>(v)];
to[static_cast<std::size_t>(edge_ptr)] = v;
next_edge[static_cast<std::size_t>(edge_ptr)] = head[static_cast<std::size_t>(u)];
head[static_cast<std::size_t>(u)] = edge_ptr;
++edge_ptr;
}
std::vector<std::uint8_t> color(static_cast<std::size_t>(N), 0U);
std::vector<int> pos(static_cast<std::size_t>(N), -1);
std::vector<std::uint8_t> in_cycle(static_cast<std::size_t>(N), 0U);
std::vector<std::vector<int>> cycles;
cycles.reserve(N / 4);
std::vector<int> path;
path.reserve(256);
for (int start = 0; start < N; ++start) {
if (color[static_cast<std::size_t>(start)] != 0U) {
continue;
}
path.clear();
int v = start;
while (color[static_cast<std::size_t>(v)] == 0U) {
color[static_cast<std::size_t>(v)] = 1U;
pos[static_cast<std::size_t>(v)] = static_cast<int>(path.size());
path.push_back(v);
v = succ[static_cast<std::size_t>(v)];
}
if (color[static_cast<std::size_t>(v)] == 1U) {
const int cycle_start = pos[static_cast<std::size_t>(v)];
std::vector<int> cyc;
cyc.reserve(path.size() - static_cast<std::size_t>(cycle_start));
for (int i = cycle_start; i < static_cast<int>(path.size()); ++i) {
const int node = path[static_cast<std::size_t>(i)];
in_cycle[static_cast<std::size_t>(node)] = 1U;
cyc.push_back(node);
}
cycles.push_back(std::move(cyc));
}
for (const int node : path) {
color[static_cast<std::size_t>(node)] = 2U;
pos[static_cast<std::size_t>(node)] = -1;
}
}
std::vector<u64> dp0(static_cast<std::size_t>(N), 0ULL);
std::vector<u64> dp1(static_cast<std::size_t>(N), 0ULL);
std::vector<std::uint8_t> done(static_cast<std::size_t>(N), 0U);
auto compute_tree = [&](const int root) {
if (done[static_cast<std::size_t>(root)] != 0U) {
return;
}
std::vector<std::pair<int, bool>> st;
st.reserve(128);
st.push_back({root, false});
while (!st.empty()) {
const auto [u, expanded] = st.back();
st.pop_back();
if (!expanded) {
st.push_back({u, true});
for (int e = head[static_cast<std::size_t>(u)]; e != -1;
e = next_edge[static_cast<std::size_t>(e)]) {
const int ch = to[static_cast<std::size_t>(e)];
if (in_cycle[static_cast<std::size_t>(ch)] != 0U) {
continue;
}
if (done[static_cast<std::size_t>(ch)] == 0U) {
st.push_back({ch, false});
}
}
continue;
}
u64 ways0 = 1ULL;
u64 ways1 = 1ULL;
for (int e = head[static_cast<std::size_t>(u)]; e != -1;
e = next_edge[static_cast<std::size_t>(e)]) {
const int ch = to[static_cast<std::size_t>(e)];
if (in_cycle[static_cast<std::size_t>(ch)] != 0U) {
continue;
}
const u64 sum = add_mod(dp0[static_cast<std::size_t>(ch)], dp1[static_cast<std::size_t>(ch)]);
ways0 = mul_mod(ways0, sum);
ways1 = mul_mod(ways1, dp0[static_cast<std::size_t>(ch)]);
}
dp0[static_cast<std::size_t>(u)] = ways0;
dp1[static_cast<std::size_t>(u)] = ways1;
done[static_cast<std::size_t>(u)] = 1U;
}
};
u64 answer = 1ULL;
for (const auto& cyc : cycles) {
const int m = static_cast<int>(cyc.size());
std::vector<u64> w0(static_cast<std::size_t>(m), 1ULL);
std::vector<u64> w1(static_cast<std::size_t>(m), 1ULL);
for (int i = 0; i < m; ++i) {
const int u = cyc[static_cast<std::size_t>(i)];
u64 a0 = 1ULL;
u64 a1 = 1ULL;
for (int e = head[static_cast<std::size_t>(u)]; e != -1;
e = next_edge[static_cast<std::size_t>(e)]) {
const int ch = to[static_cast<std::size_t>(e)];
if (in_cycle[static_cast<std::size_t>(ch)] != 0U) {
continue;
}
compute_tree(ch);
a0 = mul_mod(a0, add_mod(dp0[static_cast<std::size_t>(ch)], dp1[static_cast<std::size_t>(ch)]));
a1 = mul_mod(a1, dp0[static_cast<std::size_t>(ch)]);
}
w0[static_cast<std::size_t>(i)] = a0;
w1[static_cast<std::size_t>(i)] = a1;
}
u64 component_count = 0ULL;
if (m == 1) {
component_count = w0[0];
} else {
u64 dp_prev0 = w0[0];
u64 dp_prev1 = 0ULL;
for (int i = 1; i < m; ++i) {
const u64 ndp0 = mul_mod(add_mod(dp_prev0, dp_prev1), w0[static_cast<std::size_t>(i)]);
const u64 ndp1 = mul_mod(dp_prev0, w1[static_cast<std::size_t>(i)]);
dp_prev0 = ndp0;
dp_prev1 = ndp1;
}
component_count = add_mod(component_count, add_mod(dp_prev0, dp_prev1));
dp_prev0 = 0ULL;
dp_prev1 = w1[0];
for (int i = 1; i < m; ++i) {
const u64 ndp0 = mul_mod(add_mod(dp_prev0, dp_prev1), w0[static_cast<std::size_t>(i)]);
const u64 ndp1 = mul_mod(dp_prev0, w1[static_cast<std::size_t>(i)]);
dp_prev0 = ndp0;
dp_prev1 = ndp1;
}
component_count = add_mod(component_count, dp_prev0);
}
answer = mul_mod(answer, component_count);
}
return answer;
}
} // namespace
int main() {
assert(solve(3) == 35ULL);
assert(solve(4) == 2118ULL);
std::cout << solve(20) << '\n';
return 0;
}
Python
import sys
# Increase recursion depth for deep trees
sys.setrecursionlimit(2000000)
def solve():
n = 20
N = 1 << n
kMod = 1001001011
succ = [0] * N
for x in range(N):
b1 = x & 1
b2 = (x >> 1) & 1
b3 = (x >> 2) & 1
last = b1 & (b2 ^ b3)
succ[x] = (x >> 1) | (last << (n - 1))
head = [-1] * N
to = [0] * N
next_edge = [-1] * N
# We build the reverse graph to find trees attaching to cycles
for v in range(N):
u = succ[v]
to[v] = v
next_edge[v] = head[u]
head[u] = v
color = [0] * N
pos = [-1] * N
in_cycle = [0] * N
cycles = []
for start in range(N):
if color[start] != 0:
continue
path = []
v = start
while color[v] == 0:
color[v] = 1
pos[v] = len(path)
path.append(v)
v = succ[v]
if color[v] == 1:
cycle_start = pos[v]
cyc = []
for i in range(cycle_start, len(path)):
node = path[i]
in_cycle[node] = 1
cyc.append(node)
cycles.append(cyc)
for node in path:
color[node] = 2
pos[node] = -1
dp0 = [0] * N
dp1 = [0] * N
done = [0] * N
def compute_tree(root):
if done[root]:
return
st = [(root, False)]
while st:
u, expanded = st.pop()
if not expanded:
st.append((u, True))
e = head[u]
while e != -1:
ch = to[e]
if not in_cycle[ch]:
if not done[ch]:
st.append((ch, False))
e = next_edge[e]
else:
ways0 = 1
ways1 = 1
e = head[u]
while e != -1:
ch = to[e]
if not in_cycle[ch]:
sm = (dp0[ch] + dp1[ch]) % kMod
ways0 = (ways0 * sm) % kMod
ways1 = (ways1 * dp0[ch]) % kMod
e = next_edge[e]
dp0[u] = ways0
dp1[u] = ways1
done[u] = 1
answer = 1
for cyc in cycles:
m = len(cyc)
w0 = [1] * m
w1 = [1] * m
for i in range(m):
u = cyc[i]
a0 = 1
a1 = 1
e = head[u]
while e != -1:
ch = to[e]
if not in_cycle[ch]:
compute_tree(ch)
a0 = (a0 * (dp0[ch] + dp1[ch])) % kMod
a1 = (a1 * dp0[ch]) % kMod
e = next_edge[e]
w0[i] = a0
w1[i] = a1
component_count = 0
if m == 1:
component_count = w0[0]
else:
dp_prev0 = w0[0]
dp_prev1 = 0
for i in range(1, m):
ndp0 = ((dp_prev0 + dp_prev1) * w0[i]) % kMod
ndp1 = (dp_prev0 * w1[i]) % kMod
dp_prev0 = ndp0
dp_prev1 = ndp1
component_count = (component_count + dp_prev0 + dp_prev1) % kMod
dp_prev0 = 0
dp_prev1 = w1[0]
for i in range(1, m):
ndp0 = ((dp_prev0 + dp_prev1) * w0[i]) % kMod
ndp1 = (dp_prev0 * w1[i]) % kMod
dp_prev0 = ndp0
dp_prev1 = ndp1
component_count = (component_count + dp_prev0) % kMod
answer = (answer * component_count) % kMod
return str(answer)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
public class Euler703 {
static final long kMod = 1001001011L;
public static String solve() {
int n = 20;
int N = 1 << n;
int[] succ = new int[N];
for (int x = 0; x < N; ++x) {
int b1 = x & 1;
int b2 = (x >> 1) & 1;
int b3 = (x >> 2) & 1;
int last = b1 & (b2 ^ b3);
succ[x] = (x >> 1) | (last << (n - 1));
}
int[] head = new int[N];
Arrays.fill(head, -1);
int[] to = new int[N];
int[] nextEdge = new int[N];
for (int v = 0; v < N; ++v) {
int u = succ[v];
to[v] = v;
nextEdge[v] = head[u];
head[u] = v;
}
byte[] color = new byte[N];
int[] pos = new int[N];
Arrays.fill(pos, -1);
byte[] inCycle = new byte[N];
List<int[]> cycles = new ArrayList<>(N / 4);
int[] path = new int[256];
int pathLen = 0;
for (int start = 0; start < N; ++start) {
if (color[start] != 0) {
continue;
}
pathLen = 0;
int v = start;
while (color[v] == 0) {
color[v] = 1;
pos[v] = pathLen;
if (pathLen == path.length) {
path = Arrays.copyOf(path, path.length * 2);
}
path[pathLen++] = v;
v = succ[v];
}
if (color[v] == 1) {
int cycleStart = pos[v];
int[] cyc = new int[pathLen - cycleStart];
for (int i = cycleStart; i < pathLen; ++i) {
int node = path[i];
inCycle[node] = 1;
cyc[i - cycleStart] = node;
}
cycles.add(cyc);
}
for (int i = 0; i < pathLen; i++) {
int node = path[i];
color[node] = 2;
pos[node] = -1;
}
}
long[] dp0 = new long[N];
long[] dp1 = new long[N];
byte[] done = new byte[N];
long answer = 1;
int[] st = new int[N];
boolean[] expanded = new boolean[N];
for (int[] cyc : cycles) {
int m = cyc.length;
long[] w0 = new long[m];
long[] w1 = new long[m];
Arrays.fill(w0, 1L);
Arrays.fill(w1, 1L);
for (int i = 0; i < m; ++i) {
int startU = cyc[i];
long a0 = 1;
long a1 = 1;
for (int e = head[startU]; e != -1; e = nextEdge[e]) {
int childVec = to[e];
if (inCycle[childVec] != 0)
continue;
if (done[childVec] == 0) {
int stPtr = 0;
st[stPtr] = childVec;
expanded[stPtr] = false;
stPtr++;
while (stPtr > 0) {
stPtr--;
int u = st[stPtr];
boolean exp = expanded[stPtr];
if (!exp) {
st[stPtr] = u;
expanded[stPtr] = true;
stPtr++;
for (int ce = head[u]; ce != -1; ce = nextEdge[ce]) {
int ch = to[ce];
if (inCycle[ch] == 0 && done[ch] == 0) {
st[stPtr] = ch;
expanded[stPtr] = false;
stPtr++;
}
}
} else {
long ways0 = 1;
long ways1 = 1;
for (int ce = head[u]; ce != -1; ce = nextEdge[ce]) {
int ch = to[ce];
if (inCycle[ch] == 0) {
long sum = (dp0[ch] + dp1[ch]) % kMod;
ways0 = (ways0 * sum) % kMod;
ways1 = (ways1 * dp0[ch]) % kMod;
}
}
dp0[u] = ways0;
dp1[u] = ways1;
done[u] = 1;
}
}
}
a0 = (a0 * ((dp0[childVec] + dp1[childVec]) % kMod)) % kMod;
a1 = (a1 * dp0[childVec]) % kMod;
}
w0[i] = a0;
w1[i] = a1;
}
long componentCount = 0;
if (m == 1) {
componentCount = w0[0];
} else {
long dpPrev0 = w0[0];
long dpPrev1 = 0;
for (int i = 1; i < m; ++i) {
long ndp0 = ((dpPrev0 + dpPrev1) % kMod * w0[i]) % kMod;
long ndp1 = (dpPrev0 * w1[i]) % kMod;
dpPrev0 = ndp0;
dpPrev1 = ndp1;
}
componentCount = (componentCount + dpPrev0 + dpPrev1) % kMod;
dpPrev0 = 0;
dpPrev1 = w1[0];
for (int i = 1; i < m; ++i) {
long ndp0 = ((dpPrev0 + dpPrev1) % kMod * w0[i]) % kMod;
long ndp1 = (dpPrev0 * w1[i]) % kMod;
dpPrev0 = ndp0;
dpPrev1 = ndp1;
}
componentCount = (componentCount + dpPrev0) % kMod;
}
answer = (answer * componentCount) % kMod;
}
return Long.toString(answer);
}
public static void main(String[] args) {
System.out.println(solve());
}
}