Problem 434: Rigid Graphs
View on Project EulerProject Euler Problem 434 Solution
EulerSolve provides an optimized solution for Project Euler Problem 434, Rigid Graphs, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Consider an \(m\times n\) rectangular grid of unit squares. Each cell may be braced by a diagonal, and \(R(m,n)\) denotes the number of brace selections that make the whole framework rigid. The final target is $$S(N)=\sum_{m=1}^{N}\sum_{n=1}^{N} R(m,n)\pmod{1{,}000{,}000{,}033}.$$ The implementations do not attack rigidity directly as a geometric problem. Instead, they use the standard combinatorial model in which rigidity becomes a connectivity question in a bipartite graph. Mathematical Approach Step 1: Convert the braced grid into a bipartite graph Create one vertex for each horizontal strip of cells and one vertex for each vertical strip of cells. Thus an \(m\times n\) grid gives a bipartite graph with row vertices \(u_1,\dots,u_m\) and column vertices \(v_1,\dots,v_n\). The cell in row \(a\) and column \(b\) corresponds to the possible edge \(u_a v_b\). Choosing a set of braces is therefore equivalent to choosing a spanning subgraph of the complete bipartite graph \(K_{m,n}\). The key fact used by the solver is: $$\text{the braced framework is rigid} \iff \text{the associated bipartite graph is connected}.$$ So \(R(m,n)\) is exactly the number of connected spanning subgraphs of \(K_{m,n}\). Step 2: Count all configurations first There are \(mn\) possible row-column edges, and each edge is either present or absent....
Detailed mathematical approach
Problem Summary
Consider an \(m\times n\) rectangular grid of unit squares. Each cell may be braced by a diagonal, and \(R(m,n)\) denotes the number of brace selections that make the whole framework rigid. The final target is
$$S(N)=\sum_{m=1}^{N}\sum_{n=1}^{N} R(m,n)\pmod{1{,}000{,}000{,}033}.$$
The implementations do not attack rigidity directly as a geometric problem. Instead, they use the standard combinatorial model in which rigidity becomes a connectivity question in a bipartite graph.
Mathematical Approach
Step 1: Convert the braced grid into a bipartite graph
Create one vertex for each horizontal strip of cells and one vertex for each vertical strip of cells. Thus an \(m\times n\) grid gives a bipartite graph with row vertices \(u_1,\dots,u_m\) and column vertices \(v_1,\dots,v_n\).
The cell in row \(a\) and column \(b\) corresponds to the possible edge \(u_a v_b\). Choosing a set of braces is therefore equivalent to choosing a spanning subgraph of the complete bipartite graph \(K_{m,n}\).
The key fact used by the solver is:
$$\text{the braced framework is rigid} \iff \text{the associated bipartite graph is connected}.$$
So \(R(m,n)\) is exactly the number of connected spanning subgraphs of \(K_{m,n}\).
Step 2: Count all configurations first
There are \(mn\) possible row-column edges, and each edge is either present or absent. Therefore the total number of subgraphs is
$$2^{mn}.$$
To get \(R(m,n)\), we subtract the disconnected subgraphs.
Step 3: Classify a disconnected graph by the component of one fixed row vertex
Fix the distinguished row vertex \(u_1\). In a disconnected graph, the connected component containing \(u_1\) uses some number \(i\) of row vertices and some number \(j\) of column vertices, where
$$1\le i\le m,\qquad 0\le j\le n,\qquad (i,j)\ne(m,n).$$
For such a component:
$$\binom{m-1}{i-1}$$
chooses the other row vertices that join \(u_1\),
$$\binom{n}{j}$$
chooses the participating column vertices,
$$R(i,j)$$
counts the connected bipartite graph inside that component, and
$$2^{(m-i)(n-j)}$$
counts the arbitrary edges among the remaining rows and remaining columns. No edge is allowed between the chosen component and the remaining vertices, otherwise the component containing \(u_1\) would be larger.
Summing this contribution over all proper pairs \((i,j)\) gives the total number of disconnected graphs.
Step 4: The recurrence used by the implementations
Subtracting those disconnected cases from \(2^{mn}\) yields the split recurrence actually evaluated by the implementations, with \(M=1{,}000{,}000{,}033\):
$$R(m,n)=2^{mn}-\sum_{i=1}^{m-1}\binom{m-1}{i-1}\sum_{j=0}^{n}\binom{n}{j}R(i,j)2^{(m-i)(n-j)}-\sum_{j=0}^{n-1}\binom{n}{j}R(m,j)\pmod{M}.$$
The second sum is simply the case \(i=m\), where the factor \(2^{(m-i)(n-j)}\) becomes \(1\).
The boundary values are also important. A graph with one row vertex and no column vertices is connected, so
$$R(1,0)=1.$$
For \(m>1\), a graph with no column vertices cannot be connected, so \(R(m,0)=0\). From the recurrence one also gets \(R(1,n)=1\) for every \(n\ge 1\): with only one row vertex, all \(n\) column vertices must attach to it.
Step 5: Small worked examples
For \(R(2,2)\), start from all \(2^4=16\) subgraphs. The component containing the first row can have size \((1,0)\), \((1,1)\), or \((1,2)\), giving
$$\binom{1}{0}\left[\binom{2}{0}R(1,0)2^2+\binom{2}{1}R(1,1)2^1+\binom{2}{2}R(1,2)2^0\right]=4+4+1=9.$$
The case \(i=2\) but \(j<2\) contributes
$$\binom{2}{0}R(2,0)+\binom{2}{1}R(2,1)=0+2=2.$$
Therefore
$$R(2,2)=16-9-2=5.$$
For \(R(2,3)\), the same recurrence gives
$$\binom{1}{0}\left[\binom{3}{0}R(1,0)2^3+\binom{3}{1}R(1,1)2^2+\binom{3}{2}R(1,2)2^1+\binom{3}{3}R(1,3)2^0\right]=8+12+6+1=27,$$
and the \(i=2,\ j<3\) part equals
$$\binom{3}{0}R(2,0)+\binom{3}{1}R(2,1)+\binom{3}{2}R(2,2)=0+3+15=18.$$
Hence
$$R(2,3)=2^6-27-18=64-45=19.$$
These small values are useful sanity checks for the table.
How the Code Works
The C++, Python, and Java implementations precompute two ingredients modulo \(M\): all binomial coefficients \(\binom{a}{b}\) for \(0\le a,b\le 100\), and all powers \(2^e\) for \(0\le e\le 100^2\). After that they fill a table of \(R(m,n)\) in increasing order of \(m\) and \(n\).
For a fixed \(m\), the contribution of all terms with smaller row count \(i<m\) is accumulated first for every \(n\). Then, for each \(n\), the contribution with the same row count but smaller column count \(j<n\) is subtracted. What remains is exactly \(R(m,n)\).
The C++ implementation optionally parallelizes the accumulation over the smaller row counts, while the Python and Java implementations evaluate the same recurrence directly with ordinary nested loops. Once the table is complete, the final answer is the double sum of all \(R(m,n)\) with \(1\le m,n\le N\).
Complexity Analysis
Let the final bound be \(N\). Building Pascal's triangle and the powers of two each costs \(O(N^2)\) time. The dynamic-programming table is dominated by the triple summation over \(m\), the smaller row size \(i\), and the inner column index \(j\), which yields overall \(O(N^4)\) time. The table of values, binomial coefficients, and power table together use \(O(N^2)\) memory.
For the required case \(N=100\), this is entirely practical.
Footnotes and References
- Problem page: https://projecteuler.net/problem=434
- Bipartite graph: Wikipedia - Bipartite graph
- Connected graph: Wikipedia - Connectivity (graph theory)
- Binomial coefficient: Wikipedia - Binomial coefficient
Problem 434 source code
C++
#include <algorithm>
#include <atomic>
#include <chrono>
#include <cstdint>
#include <iostream>
#include <string>
#include <thread>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
constexpr u64 kMod = 1'000'000'033ULL;
constexpr int kDefaultN = 100;
u64 add_mod(u64 a, u64 b) {
a += b;
if (a >= kMod) {
a -= kMod;
}
return a;
}
u64 sub_mod(u64 a, u64 b) {
return (a >= b) ? (a - b) : (a + kMod - b);
}
u64 mul_mod(u64 a, u64 b) {
return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % kMod);
}
unsigned choose_thread_count(bool allow_multithreading,
unsigned requested_threads,
std::size_t workload) {
if (!allow_multithreading || workload < 300'000) {
return 1;
}
unsigned threads = requested_threads;
if (threads == 0) {
threads = std::thread::hardware_concurrency();
if (threads == 0) {
threads = 1;
}
}
threads = std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
return threads;
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
unsigned& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
unsigned parsed = 0;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10 + static_cast<unsigned>(c - '0');
}
value = parsed;
return true;
}
bool parse_int_after_prefix(const std::string& arg, const char* prefix, int& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
int parsed = 0;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10 + (c - '0');
}
value = parsed;
return true;
}
u64 brute_connected_bipartite(int m, int n) {
const int edges = m * n;
if (edges < 0 || edges > 20) {
return 0;
}
const int vertices = m + n;
const u64 masks = 1ULL << edges;
std::vector<int> queue(vertices);
std::vector<char> seen(vertices, 0);
u64 connected_count = 0;
for (u64 mask = 0; mask < masks; ++mask) {
std::fill(seen.begin(), seen.end(), 0);
int head = 0;
int tail = 0;
seen[0] = 1;
queue[tail++] = 0;
while (head < tail) {
const int u = queue[head++];
if (u < m) {
const int row = u;
for (int col = 0; col < n; ++col) {
const int edge_index = row * n + col;
if (((mask >> edge_index) & 1ULL) == 0ULL) {
continue;
}
const int v = m + col;
if (!seen[v]) {
seen[v] = 1;
queue[tail++] = v;
}
}
} else {
const int col = u - m;
for (int row = 0; row < m; ++row) {
const int edge_index = row * n + col;
if (((mask >> edge_index) & 1ULL) == 0ULL) {
continue;
}
const int v = row;
if (!seen[v]) {
seen[v] = 1;
queue[tail++] = v;
}
}
}
}
bool connected = true;
for (const char s : seen) {
if (!s) {
connected = false;
break;
}
}
if (connected) {
++connected_count;
}
}
return connected_count;
}
struct Options {
int n = kDefaultN;
bool allow_multithreading = true;
unsigned requested_threads = 0;
};
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--single-thread") {
options.allow_multithreading = false;
continue;
}
unsigned threads = 0;
if (parse_unsigned_after_prefix(arg, "--threads=", threads)) {
options.requested_threads = threads;
continue;
}
int n = 0;
if (parse_int_after_prefix(arg, "--n=", n)) {
if (n < 1) {
std::cerr << "Invalid --n value: " << n << '\n';
return false;
}
options.n = n;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return true;
}
class Euler434Solver {
public:
Euler434Solver(int n, bool allow_multithreading, unsigned requested_threads)
: n_(n),
allow_multithreading_(allow_multithreading),
requested_threads_(requested_threads),
binom_(static_cast<std::size_t>(n_) + 1,
std::vector<u64>(static_cast<std::size_t>(n_) + 1, 0)),
pow2_(static_cast<std::size_t>(n_ * n_) + 1, 0),
r_(static_cast<std::size_t>(n_) + 1,
std::vector<u64>(static_cast<std::size_t>(n_) + 1, 0)) {
build_binomials();
build_powers_of_two();
}
void solve() {
if (n_ >= 1) {
r_[1][0] = 1; // Connected component with one fixed row vertex and no column.
}
for (int m = 1; m <= n_; ++m) {
const std::vector<u64> base_prev = compute_base_from_previous_rows(m);
for (int n = 1; n <= n_; ++n) {
u64 same_row = 0;
for (int j = 0; j < n; ++j) {
same_row = add_mod(same_row,
mul_mod(binom_[static_cast<std::size_t>(n)]
[static_cast<std::size_t>(j)],
r_[static_cast<std::size_t>(m)]
[static_cast<std::size_t>(j)]));
}
u64 value = pow2_[static_cast<std::size_t>(m * n)];
value = sub_mod(value, base_prev[static_cast<std::size_t>(n)]);
value = sub_mod(value, same_row);
r_[static_cast<std::size_t>(m)][static_cast<std::size_t>(n)] = value;
}
}
}
u64 r_value(int m, int n) const {
if (m < 0 || m > n_ || n < 0 || n > n_) {
return 0;
}
return r_[static_cast<std::size_t>(m)][static_cast<std::size_t>(n)];
}
u64 s_value(int n) const {
if (n < 1) {
return 0;
}
n = std::min(n, n_);
u64 total = 0;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
total = add_mod(total,
r_[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)]);
}
}
return total;
}
bool run_checkpoints() const {
struct Checkpoint {
int m;
int n;
u64 expected;
const char* label;
};
constexpr Checkpoint checks[] = {
{2, 3, 19, "R(2,3)"},
{5, 5, 23'679'901ULL, "R(5,5)"},
};
for (const Checkpoint& checkpoint : checks) {
if (checkpoint.m > n_ || checkpoint.n > n_) {
continue;
}
const u64 got = r_value(checkpoint.m, checkpoint.n);
if (got != checkpoint.expected) {
std::cerr << "Checkpoint failed for " << checkpoint.label
<< ": expected " << checkpoint.expected
<< ", got " << got << '\n';
return false;
}
std::cout << checkpoint.label << " = " << got << " (ok)" << '\n';
}
if (n_ >= 5) {
constexpr u64 expected_s5 = 25'021'721ULL;
const u64 got_s5 = s_value(5);
if (got_s5 != expected_s5) {
std::cerr << "Checkpoint failed for S(5): expected " << expected_s5
<< ", got " << got_s5 << '\n';
return false;
}
std::cout << "S(5) = " << got_s5 << " (ok)" << '\n';
}
struct BruteCheckpoint {
int m;
int n;
};
constexpr BruteCheckpoint brute_checks[] = {
{2, 2},
{2, 3},
{3, 3},
};
for (const BruteCheckpoint& checkpoint : brute_checks) {
if (checkpoint.m > n_ || checkpoint.n > n_) {
continue;
}
const u64 brute = brute_connected_bipartite(checkpoint.m, checkpoint.n);
const u64 got = r_value(checkpoint.m, checkpoint.n);
if (got != brute) {
std::cerr << "Brute-force checkpoint failed for R(" << checkpoint.m
<< ',' << checkpoint.n << "): expected " << brute
<< ", got " << got << '\n';
return false;
}
std::cout << "Brute R(" << checkpoint.m << ',' << checkpoint.n << ") = " << got
<< " (ok)" << '\n';
}
if (n_ >= 5) {
const u64 lhs = r_value(3, 5);
const u64 rhs = r_value(5, 3);
if (lhs != rhs) {
std::cerr << "Symmetry checkpoint failed: R(3,5)=" << lhs
<< ", R(5,3)=" << rhs << '\n';
return false;
}
std::cout << "Symmetry R(3,5)=R(5,3) = " << lhs << " (ok)" << '\n';
}
return true;
}
private:
int n_ = 0;
bool allow_multithreading_ = true;
unsigned requested_threads_ = 0;
std::vector<std::vector<u64>> binom_;
std::vector<u64> pow2_;
std::vector<std::vector<u64>> r_;
void build_binomials() {
binom_[0][0] = 1;
for (int n = 1; n <= n_; ++n) {
binom_[static_cast<std::size_t>(n)][0] = 1;
binom_[static_cast<std::size_t>(n)][static_cast<std::size_t>(n)] = 1;
for (int k = 1; k < n; ++k) {
binom_[static_cast<std::size_t>(n)][static_cast<std::size_t>(k)] =
add_mod(binom_[static_cast<std::size_t>(n - 1)]
[static_cast<std::size_t>(k - 1)],
binom_[static_cast<std::size_t>(n - 1)]
[static_cast<std::size_t>(k)]);
}
}
}
void build_powers_of_two() {
pow2_[0] = 1;
for (int e = 1; e <= n_ * n_; ++e) {
pow2_[static_cast<std::size_t>(e)] =
add_mod(pow2_[static_cast<std::size_t>(e - 1)],
pow2_[static_cast<std::size_t>(e - 1)]);
}
}
void accumulate_previous_row_contribution(int m,
int i,
std::vector<u64>& out) const {
const u64 choose_rows = binom_[static_cast<std::size_t>(m - 1)]
[static_cast<std::size_t>(i - 1)];
const int delta = m - i;
for (int n = 1; n <= n_; ++n) {
u64 inner = 0;
for (int j = 0; j <= n; ++j) {
const u64 choose_cols =
binom_[static_cast<std::size_t>(n)][static_cast<std::size_t>(j)];
const u64 rij = r_[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)];
const u64 free_part =
pow2_[static_cast<std::size_t>(delta * (n - j))];
const u64 term = mul_mod(mul_mod(choose_cols, rij), free_part);
inner = add_mod(inner, term);
}
out[static_cast<std::size_t>(n)] =
add_mod(out[static_cast<std::size_t>(n)], mul_mod(choose_rows, inner));
}
}
std::vector<u64> compute_base_from_previous_rows(int m) const {
std::vector<u64> base(static_cast<std::size_t>(n_) + 1, 0);
if (m <= 1) {
return base;
}
const std::size_t workload =
static_cast<std::size_t>(m - 1) * static_cast<std::size_t>(n_) *
static_cast<std::size_t>(n_ + 1) / 2;
const unsigned threads =
choose_thread_count(allow_multithreading_, requested_threads_, workload);
if (threads == 1) {
for (int i = 1; i < m; ++i) {
accumulate_previous_row_contribution(m, i, base);
}
return base;
}
std::atomic<int> next_i{1};
std::vector<std::vector<u64>> partial(
threads, std::vector<u64>(static_cast<std::size_t>(n_) + 1, 0));
std::vector<std::thread> workers;
workers.reserve(threads);
for (unsigned tid = 0; tid < threads; ++tid) {
workers.emplace_back([&, tid]() {
while (true) {
const int i = next_i.fetch_add(1, std::memory_order_relaxed);
if (i >= m) {
break;
}
accumulate_previous_row_contribution(m, i, partial[tid]);
}
});
}
for (std::thread& worker : workers) {
worker.join();
}
for (unsigned tid = 0; tid < threads; ++tid) {
for (int n = 1; n <= n_; ++n) {
base[static_cast<std::size_t>(n)] =
add_mod(base[static_cast<std::size_t>(n)],
partial[tid][static_cast<std::size_t>(n)]);
}
}
return base;
}
};
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
std::cerr << "Usage: ./Euler434 [--n=100] [--single-thread] [--threads=K]" << '\n';
return 1;
}
const auto t0 = std::chrono::steady_clock::now();
Euler434Solver solver(options.n,
options.allow_multithreading,
options.requested_threads);
solver.solve();
if (!solver.run_checkpoints()) {
return 1;
}
const u64 answer = solver.s_value(options.n);
const auto t1 = std::chrono::steady_clock::now();
const std::chrono::duration<double> elapsed = t1 - t0;
std::cout << "S(" << options.n << ") mod " << kMod << " = " << answer << '\n';
std::cout << "Elapsed: " << elapsed.count() << " s" << '\n';
return 0;
}
Python
def solve():
MOD = 1_000_000_033
N = 100
def add_mod(a, b): return (a + b) % MOD
def sub_mod(a, b): return (a - b) % MOD
def mul_mod(a, b): return a * b % MOD
# Build binomials
binom = [[0]*(N+1) for _ in range(N+1)]
binom[0][0] = 1
for n in range(1, N+1):
binom[n][0] = 1
binom[n][n] = 1
for k in range(1, n):
binom[n][k] = add_mod(binom[n-1][k-1], binom[n-1][k])
# Build powers of 2
pow2 = [0] * (N*N + 1)
pow2[0] = 1
for e in range(1, N*N+1):
pow2[e] = add_mod(pow2[e-1], pow2[e-1])
r = [[0]*(N+1) for _ in range(N+1)]
r[1][0] = 1
for m in range(1, N+1):
# Compute base from previous rows
base = [0] * (N + 1)
for i in range(1, m):
choose_rows = binom[m-1][i-1]
delta = m - i
for n in range(1, N+1):
inner = 0
for j in range(n+1):
term = mul_mod(mul_mod(binom[n][j], r[i][j]), pow2[delta*(n-j)])
inner = add_mod(inner, term)
base[n] = add_mod(base[n], mul_mod(choose_rows, inner))
for n in range(1, N+1):
same_row = 0
for j in range(n):
same_row = add_mod(same_row, mul_mod(binom[n][j], r[m][j]))
value = pow2[m * n]
value = sub_mod(value, base[n])
value = sub_mod(value, same_row)
r[m][n] = value
total = 0
for i in range(1, N+1):
for j in range(1, N+1):
total = add_mod(total, r[i][j])
return str(total)
if __name__ == '__main__':
print(solve())
Java
public class Euler434 {
private static final long MOD = 1000000033L;
private static final int N_VAL = 100;
public static String solve() {
long[][] binom = new long[N_VAL + 1][N_VAL + 1];
binom[0][0] = 1;
for (int i = 1; i <= N_VAL; i++) {
binom[i][0] = 1;
binom[i][i] = 1;
for (int j = 1; j < i; j++) {
binom[i][j] = (binom[i - 1][j - 1] + binom[i - 1][j]) % MOD;
}
}
long[] pow2 = new long[N_VAL * N_VAL + 1];
pow2[0] = 1;
for (int i = 1; i <= N_VAL * N_VAL; i++) {
pow2[i] = (pow2[i - 1] * 2) % MOD;
}
long[][] r = new long[N_VAL + 1][N_VAL + 1];
r[1][0] = 1;
for (int m = 1; m <= N_VAL; m++) {
long[] base = new long[N_VAL + 1];
if (m > 1) {
for (int i = 1; i < m; i++) {
long chooseRows = binom[m - 1][i - 1];
int delta = m - i;
for (int n = 1; n <= N_VAL; n++) {
long inner = 0;
for (int j = 0; j <= n; j++) {
long term = (binom[n][j] * r[i][j]) % MOD;
term = (term * pow2[delta * (n - j)]) % MOD;
inner = (inner + term) % MOD;
}
base[n] = (base[n] + chooseRows * inner) % MOD;
}
}
}
for (int n = 1; n <= N_VAL; n++) {
long sameRow = 0;
for (int j = 0; j < n; j++) {
sameRow = (sameRow + binom[n][j] * r[m][j]) % MOD;
}
long val = pow2[m * n];
val = (val - base[n]) % MOD;
if (val < 0)
val += MOD;
val = (val - sameRow) % MOD;
if (val < 0)
val += MOD;
r[m][n] = val;
}
}
long total = 0;
for (int i = 1; i <= N_VAL; i++) {
for (int j = 1; j <= N_VAL; j++) {
total = (total + r[i][j]) % MOD;
}
}
return String.valueOf(total);
}
public static void main(String[] args) {
System.out.println(solve());
}
}