Problem 701: Random Connected Area
View on Project EulerProject Euler Problem 701 Solution
EulerSolve provides an optimized solution for Project Euler Problem 701, Random Connected Area, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary We consider a random \(W\times H\) binary grid, with the Project Euler instance using \(W=H=7\). Each cell is occupied independently with probability \(1/2\), and occupied cells are connected by 4-neighbor adjacency. If \(L(B)\) denotes the size of the largest occupied connected component of a board \(B\), the task is to compute the exact expectation \(\mathbb E[L(B)]\). A brute-force search would inspect all \(2^{WH}\) boards, which is impossible for \(7\times 7\). The successful idea is to scan the grid one column at a time and keep only the information that can still influence future growth. Mathematical Approach Step 1: Rewrite the expectation as an exact counting problem Let $$\Omega_{W,H}=\{0,1\}^{W\times H}.$$ For each board \(B\in\Omega_{W,H}\), let \(\mathcal C(B)\) be the set of occupied 4-connected components of \(B\), and define $$L(B)=\max_{C\in\mathcal C(B)} \lvert C\rvert,$$ with the convention \(L(B)=0\) when the board has no occupied cells. Since every board is equally likely, $$\mathbb E[L]=\frac{1}{2^{WH}}\sum_{B\in\Omega_{W,H}} L(B).$$ So the whole problem is to sum the largest-component size over all boards exactly, then divide by \(2^{WH}\). Step 2: Keep only the frontier information that matters Suppose the first \(c\) columns have already been processed....
Detailed mathematical approach
Problem Summary
We consider a random \(W\times H\) binary grid, with the Project Euler instance using \(W=H=7\). Each cell is occupied independently with probability \(1/2\), and occupied cells are connected by 4-neighbor adjacency. If \(L(B)\) denotes the size of the largest occupied connected component of a board \(B\), the task is to compute the exact expectation \(\mathbb E[L(B)]\).
A brute-force search would inspect all \(2^{WH}\) boards, which is impossible for \(7\times 7\). The successful idea is to scan the grid one column at a time and keep only the information that can still influence future growth.
Mathematical Approach
Step 1: Rewrite the expectation as an exact counting problem
Let
$$\Omega_{W,H}=\{0,1\}^{W\times H}.$$
For each board \(B\in\Omega_{W,H}\), let \(\mathcal C(B)\) be the set of occupied 4-connected components of \(B\), and define
$$L(B)=\max_{C\in\mathcal C(B)} \lvert C\rvert,$$
with the convention \(L(B)=0\) when the board has no occupied cells. Since every board is equally likely,
$$\mathbb E[L]=\frac{1}{2^{WH}}\sum_{B\in\Omega_{W,H}} L(B).$$
So the whole problem is to sum the largest-component size over all boards exactly, then divide by \(2^{WH}\).
Step 2: Keep only the frontier information that matters
Suppose the first \(c\) columns have already been processed. Any connected component that does not touch column \(c\) can never grow again, because all future cells lie strictly to the right. Therefore closed components only matter through a single number: the largest closed size seen so far.
Now look at column \(c\). Its occupied cells form vertical runs such as \(11100 11\mapsto (111)\) and \((11)\). Let these runs be \(R_1,\dots,R_t\). Two different runs in the same column may already belong to the same global component because earlier columns can connect them around the gaps. Hence the frontier is not described only by the mask of occupied cells; we must also know which runs belong together.
The state can therefore be written as
$$\sigma=(b,\Pi,a_1,\dots,a_r),$$
where:
Here \(b\) denotes the largest closed component size seen so far.
\(\Pi\) is a partition of the runs \(R_1,\dots,R_t\), and each block of \(\Pi\) corresponds to one open component touching the frontier,
For \(1\le j\le r\), \(a_j\) denotes the total size of the \(j\)-th open component.
This is sufficient because future columns can interact with the processed region only through those open components.
Step 3: Append one new column
Let the next column be the bitmask \(m\in\{0,1\}^W\). Its occupied cells also decompose into vertical runs, say \(S_1,\dots,S_u\).
An old open component touches the new column exactly at rows where both adjacent columns contain an occupied cell. Once one cell of a run \(S_k\) is touched, the whole run belongs to the same new connected piece, because cells inside one vertical run are already connected within the new column.
This produces three possibilities:
First, an old open component may touch nothing in the new column. Then it becomes closed forever and can only update the best closed value.
Second, one old component may touch one or more runs in the new column. Then it survives and its size grows by the number of newly occupied cells attached to it.
Third, several old components may touch the same run, or runs that become linked inside the new column. Then those components merge into one new open component.
If \(\mathcal D\) is the set of old open components that disappear, then
$$b'=\max\!\Bigl(b,\max_{j\in\mathcal D} a_j\Bigr).$$
For each new open component \(T\), let \(\mathcal A(T)\) be the set of old components that feed into it. Its new size is
$$a'(T)=\lvert T\rvert+\sum_{j\in\mathcal A(T)} a_j,$$
where \(\lvert T\rvert\) counts only the occupied cells in the newly appended column.
Step 4: Why partitions of runs are enough
Inside a single column, connectivity is completely determined by vertical runs. Cells in the same run are automatically connected, while cells in different runs are separated inside that column. The only missing information is whether two distinct runs are already linked through earlier columns. That is exactly what the partition \(\Pi\) records.
This compression is powerful because a width-7 column has at most four occupied runs, for example \(1010101\). So for each column mask there are only finitely many possible partitions of those runs, and all of them can be precomputed once. The dynamic program then works with canonical frontier patterns instead of full partial boards.
Step 5: Worked example
Take width \(W=3\). Suppose the first processed column is \(101\). Then there are two runs:
$$R_1=\{1\},\qquad R_2=\{3\}.$$
No earlier column exists, so they are separate open components of sizes \(1\) and \(1\). The state is
$$\sigma=\bigl(0,\{\{R_1\},\{R_2\}\},1,1\bigr).$$
Now append the next column \(111\). In the new column there is one run
$$S_1=\{1,2,3\}.$$
The old component at row \(1\) touches \(S_1\), and the old component at row \(3\) also touches \(S_1\). Because those contacts lie in the same vertical run, the whole run becomes one connected bridge, so the two old components merge. The new open component has size
$$1+1+3=5.$$
The new state is therefore
$$\sigma'=\bigl(0,\{\{S_1\}\},5\bigr).$$
If the following column were \(000\), that size-5 component would lose frontier contact, close permanently, and update the best closed value to \(5\).
Step 6: Recover the expectation from state counts
Let \(w_c(\sigma)\) be the number of \(c\)-column prefixes that produce state \(\sigma\). The transition defined above gives
$$w_{c+1}(\sigma')=\sum_{\sigma}\sum_{\substack{m\in\{0,1\}^W\\ T_m(\sigma)=\sigma'}} w_c(\sigma),$$
where \(T_m\) is the deterministic update induced by the next column mask \(m\). The dynamic program stores these weights as exact integers, not floating-point probabilities.
After all \(H\) columns have been processed, every terminal state \(\sigma=(b,\Pi,a_1,\dots,a_r)\) contributes
$$\operatorname{score}(\sigma)=\max(b,a_1,\dots,a_r)$$
to the largest connected area of each board counted by that state. Hence
$$\mathbb E[L]=\frac{1}{2^{WH}}\sum_{\sigma\in\mathcal T} w_H(\sigma)\operatorname{score}(\sigma),$$
where \(\mathcal T\) is the set of terminal states. This is exactly the quantity evaluated by the implementation.
How the Code Works
The implementation first enumerates every possible column mask, decomposes each mask into vertical runs, and precomputes the canonical partitions of those runs. It also precomputes how an old frontier pattern interacts with a new column: which old components survive, which ones merge, which ones disappear, and how many new occupied cells are absorbed into each surviving component.
The dynamic program then advances one column at a time. A state is keyed by the frontier partition together with the accumulated open-component sizes and the current best closed size. When a new column is appended, the implementation applies the precomputed connectivity rule, creates the resulting state, and merges identical target states by adding their integer multiplicities.
Because every column mask is equally likely, the C++, Python, and Java implementations count boards rather than carrying fractions at each step. Only at the very end do they divide by \(2^{WH}\). The Python and Java implementations use the same exact underlying solver as the C++ implementation, so all three languages follow identical transitions and produce the same final value.
The implementation also uses horizontal reflection on the penultimate layer: mirror-image masks contribute equally to the final expectation, so non-symmetric representatives can be counted once with multiplicity \(2\). This reduces the amount of final work without changing the exact result.
Complexity Analysis
Let \(P\) be the total number of canonical frontier partitions across all column masks, and let \(S_c\) be the number of reachable compressed states after \(c\) columns. Precomputing mask-to-mask contact information costs \(O(4^W\cdot W)\), and building the partition-transition tables costs \(O(P\cdot 2^W)\).
The main dynamic program processes each reachable state against each possible next-column mask, so its time complexity is
$$O\!\left(\sum_{c=0}^{H-1} 2^W S_c\right).$$
The memory usage is
$$O(S_c+P\cdot 2^W)$$
for the active DP layers and the precomputed transition data. For the \(7\times 7\) instance this is practical, while direct enumeration of all \(2^{49}\) boards is not.
Footnotes and References
- Problem page: https://projecteuler.net/problem=701
- Connected components in graph theory: Wikipedia — Connected component
- Transfer-matrix method: Wikipedia — Transfer-matrix method
- Bell numbers and set partitions: Wikipedia — Bell number
- Dynamic programming: Wikipedia — Dynamic programming
Problem 701 source code
C++
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <assert.h>
#include <cstdint>
#include <pthread.h>
#include <unistd.h>
#include <vector>
#define DEFAULT_WIDTH 7
#define DEFAULT_HEIGHT 7
#define MAX_CONNECTED 5
#define USE_128BIT_WEIGHT
#if defined(USE_128BIT_WEIGHT)
using ProblemWeight = __int128 ;
#define MAX_COLUMN_LIMIT 9
#else
using ProblemWeight = uint64_t ;
#define MAX_COLUMN_LIMIT 7
#endif
typedef uint32_t IndexValue ;
typedef uint8_t SurfaceValue ;
typedef struct CombinationIndex {
IndexValue N ;
IndexValue k ;
IndexValue *index ;
} CombinationIndex ;
CombinationIndex * combinationIndexAlloc(IndexValue N, IndexValue k) ;
CombinationIndex * combinationIndexDestroy(CombinationIndex * NS) ;
IndexValue combinationIndexSize(const CombinationIndex *NS,int ks) ;
IndexValue combinationIndexRank(const CombinationIndex *NS,int ks,SurfaceValue *sum) ;
typedef struct PartitionLayout {
int id ;
int source_mask ;
int set_count ;
uint16_t part_masks[MAX_CONNECTED] ;
} PartitionLayout ;
typedef struct PartitionTransition {
int next_partition ;
uint8_t source_index[MAX_CONNECTED] ;
uint8_t added_surface[MAX_CONNECTED] ;
} PartitionTransition ;
typedef struct DpState {
uint16_t partition_id ;
SurfaceValue surface[MAX_CONNECTED+1] ;
ProblemWeight prob ;
} DpState ;
typedef struct GeneratedState {
IndexValue hash_key ;
uint16_t partition_id ;
SurfaceValue surface[MAX_CONNECTED+1] ;
ProblemWeight prob ;
} GeneratedState ;
typedef struct BuildRangeTask {
const DpState *states ;
int state_start ;
int state_end ;
const PartitionLayout *partitions ;
const PartitionTransition *transitions ;
const CombinationIndex *hash_basis ;
const IndexValue *partition_base_index ;
int mask_count ;
int added_mask ;
int duplicate_factor ;
int hash_table_size ;
std::vector<GeneratedState> locals ;
} BuildRangeTask ;
typedef struct LayerSumTask {
const DpState *states ;
const PartitionLayout *partitions ;
int state_start ;
int state_end ;
ProblemWeight partial ;
} LayerSumTask ;
typedef struct SymmetrySumTask {
const DpState *states ;
const PartitionLayout *partitions ;
const PartitionTransition *transitions ;
int mask_count ;
int state_start ;
int state_end ;
ProblemWeight partial ;
} SymmetrySumTask ;
typedef struct PrecomputeData {
int mask_count ;
int partition_count ;
int *symmetry ;
PartitionLayout *partitions ;
int *first_partition_by_mask ;
PartitionTransition *transitions ;
} PrecomputeData ;
void buildPrecompute(PrecomputeData *solver_data, int width) ;
static int get_thread_count() {
long n = sysconf(_SC_NPROCESSORS_ONLN);
if(n < 1) return 1;
return (int)n;
}
static void *buildTransitionWorker(void *arg) {
BuildRangeTask *task = (BuildRangeTask *)arg;
if(task->state_start >= task->state_end) {
return nullptr;
}
task->locals.clear();
task->locals.reserve((size_t)(task->state_end - task->state_start));
std::vector<IndexValue> localHash(task->hash_table_size, 0);
ProblemWeight addMul = task->duplicate_factor ? 2 : 1;
for(int state_index = task->state_start; state_index < task->state_end; ++state_index) {
const DpState *source_state = task->states + state_index;
int source_partition = source_state->partition_id ;
const PartitionTransition *trf = task->transitions + source_partition * task->mask_count + task->added_mask ;
int target_partition = trf->next_partition ;
SurfaceValue surf[MAX_CONNECTED+1] ;
surf[0] = source_state->surface[0] ;
for(int segment = 0; segment < task->partitions[target_partition].set_count; segment++) {
surf[segment+1] = trf->added_surface[segment] ;
}
for(int segment = 0; segment < task->partitions[source_partition].set_count; segment++) {
if(trf->source_index[segment]) {
surf[trf->source_index[segment]] += source_state->surface[segment+1];
} else {
if(source_state->surface[segment+1] > surf[0]) surf[0] = source_state->surface[segment+1] ;
}
}
IndexValue numHash = task->partition_base_index[target_partition] + combinationIndexRank(task->hash_basis, task->partitions[target_partition].set_count+1, surf) ;
IndexValue loc = localHash[numHash] ;
ProblemWeight add = addMul * source_state->prob ;
if(loc == 0) {
GeneratedState state ;
state.hash_key = numHash ;
state.partition_id = (uint16_t)target_partition ;
memcpy(state.surface, surf, sizeof(state.surface)) ;
state.prob = add ;
task->locals.push_back(state) ;
localHash[numHash] = (IndexValue) task->locals.size() ;
} else {
task->locals[loc - 1].prob += add ;
}
}
return nullptr;
}
static void *accumulateLayerWorker(void *arg) {
LayerSumTask *task = (LayerSumTask *)arg;
ProblemWeight total = 0 ;
for(int state_index = task->state_start; state_index < task->state_end; state_index++) {
const DpState *new_state = task->states + state_index ;
int maxS = new_state->surface[0];
const PartitionLayout *pp = task->partitions + new_state->partition_id ;
for(int segment = 0; segment < pp->set_count; segment++) {
if(new_state->surface[segment+1] > maxS) maxS = new_state->surface[segment+1];
}
total += (ProblemWeight)maxS * new_state->prob ;
}
task->partial = total;
return nullptr;
}
static void *accumulateSymmetryWorker(void *arg) {
SymmetrySumTask *task = (SymmetrySumTask *)arg;
ProblemWeight total = 0 ;
for(int state_index = task->state_start; state_index < task->state_end; state_index++) {
const DpState *source_state = task->states + state_index ;
int source_partition = source_state->partition_id ;
ProblemWeight state_prob = source_state->prob ;
for(int column = 0; column < task->mask_count; column++) {
const PartitionTransition *trf = task->transitions + source_partition*task->mask_count + column ;
int target_partition = trf->next_partition ;
SurfaceValue surf[MAX_CONNECTED] ;
int maxS = source_state->surface[0] ;
for(int segment = 0; segment < task->partitions[target_partition].set_count; segment++) {
surf[segment] = trf->added_surface[segment] ;
}
for(int segment = 0; segment < task->partitions[source_partition].set_count; segment++) {
if(trf->source_index[segment]) {
surf[trf->source_index[segment]-1] += source_state->surface[segment+1];
} else {
if(source_state->surface[segment+1] > maxS ) maxS = source_state->surface[segment+1] ;
}
}
for(int segment = 0; segment < task->partitions[target_partition].set_count; segment++) {
if(surf[segment] > maxS) maxS = surf[segment] ;
}
total += (ProblemWeight)maxS * state_prob;
}
}
task->partial = total;
return nullptr;
}
int main(int argc, const char * argv[]) {
PrecomputeData solver{};
int width = DEFAULT_WIDTH ;
int maxNbColumn = DEFAULT_HEIGHT ;
if(argc > 1){
width = atoi(argv[1]) ;
if(width > MAX_COLUMN_LIMIT) width = MAX_COLUMN_LIMIT ;
if(argc > 2) {
maxNbColumn = atoi(argv[2]);
if(maxNbColumn > MAX_COLUMN_LIMIT) maxNbColumn = MAX_COLUMN_LIMIT ;
if(maxNbColumn < width) {
int tmp = width ;
width = maxNbColumn ;
maxNbColumn = tmp ;
}
}
}
buildPrecompute(&solver, width) ;
int maskCount = solver.mask_count ;
const PartitionLayout *partitions = solver.partitions ;
const PartitionTransition *transitions = solver.transitions;
DpState *states = (DpState*) malloc(300000000 * sizeof(states[0])) ;
int stateCount = 1;
IndexValue *partitionBaseIndex = (IndexValue*) calloc(solver.partition_count, sizeof(partitionBaseIndex[0])) ;
int firstStateInPrevLayer = 0 ;
memset(states, 0, sizeof(states[0]));
states[0].prob = 1 ;
states[0].partition_id = 0 ;
int firstStateInCurrentLayer = 1;
double result ;
for(int columns = 1; columns < maxNbColumn; columns++) {
int totalWeight = width * columns ;
CombinationIndex *hashIndex = combinationIndexAlloc(totalWeight, MAX_CONNECTED+1) ;
ProblemWeight sumCurrent = 0 , sumSymmetric = 0 ;
for(int addedMask = 0; addedMask < maskCount; addedMask++) {
int symmetryMultiplier = 0 ;
if(columns == maxNbColumn - 1) {
if(solver.symmetry[addedMask] < 0) continue ;
symmetryMultiplier = solver.symmetry[addedMask] ? 1 : 0 ;
}
int hashSize = 0 ;
int firstStateForMask = stateCount ;
for(int partitionId = solver.first_partition_by_mask[addedMask]; partitionId < solver.first_partition_by_mask[addedMask + 1]; partitionId++) {
partitionBaseIndex[partitionId] = hashSize ;
hashSize += combinationIndexSize(hashIndex, partitions[partitionId].set_count + 1) ;
}
std::vector<pthread_t> buildThreads;
std::vector<BuildRangeTask> buildTasks;
int activeStateCount = firstStateInCurrentLayer - firstStateInPrevLayer ;
int buildThreadCount = get_thread_count();
if(activeStateCount < 4096) buildThreadCount = 1 ;
if(buildThreadCount > activeStateCount) buildThreadCount = activeStateCount;
if(buildThreadCount <= 0) buildThreadCount = 1;
buildThreads.resize(buildThreadCount);
buildTasks.resize(buildThreadCount);
int block = (activeStateCount + buildThreadCount -1) / buildThreadCount ;
for(int thread = 0; thread < buildThreadCount; thread++) {
int start = firstStateInPrevLayer + thread * block ;
int end = start + block ;
if(end > firstStateInCurrentLayer) end = firstStateInCurrentLayer;
BuildRangeTask &task = buildTasks[thread];
task.states = states;
task.state_start = start;
task.state_end = end;
task.partitions = partitions;
task.transitions = transitions;
task.hash_basis = hashIndex;
task.partition_base_index = partitionBaseIndex;
task.mask_count = maskCount;
task.added_mask = addedMask;
task.duplicate_factor = symmetryMultiplier;
task.hash_table_size = hashSize;
pthread_create(&buildThreads[thread], NULL, buildTransitionWorker, &task);
}
for(int thread = 0; thread < buildThreadCount; thread++) pthread_join(buildThreads[thread], NULL);
std::vector<IndexValue> stateSlot(hashSize, 0);
for(int thread = 0; thread < buildThreadCount; thread++) {
const std::vector<GeneratedState> &local = buildTasks[thread].locals;
for(size_t i = 0; i < local.size(); ++i) {
const GeneratedState &generated = local[i];
IndexValue existing = stateSlot[generated.hash_key];
if(existing == 0) {
stateSlot[generated.hash_key] = (IndexValue)(stateCount + 1);
DpState *nextState = states + stateCount++;
nextState->partition_id = generated.partition_id ;
nextState->prob = generated.prob ;
memcpy(nextState->surface, generated.surface, sizeof(nextState->surface));
} else {
states[existing - 1].prob += generated.prob ;
}
}
}
int statesInSlice = stateCount - firstStateForMask ;
if(statesInSlice > 0) {
int sumThreads = get_thread_count();
if(statesInSlice < 2048) sumThreads = 1 ;
if(sumThreads > statesInSlice) sumThreads = statesInSlice;
if(sumThreads <= 0) sumThreads = 1;
std::vector<pthread_t> sumWorkers(sumThreads);
std::vector<LayerSumTask> sumTasks(sumThreads);
int per = (statesInSlice + sumThreads -1) / sumThreads ;
for(int thread = 0; thread < sumThreads; thread++) {
int start = firstStateForMask + thread * per ;
int end = start + per ;
if(end > stateCount) end = stateCount;
LayerSumTask &task = sumTasks[thread];
task.states = states;
task.partitions = partitions;
task.state_start = start;
task.state_end = end;
task.partial = 0;
pthread_create(&sumWorkers[thread], NULL, accumulateLayerWorker, &task);
}
for(int thread = 0; thread < sumThreads; thread++) {
pthread_join(sumWorkers[thread], NULL);
sumCurrent += sumTasks[thread].partial ;
}
}
if(columns == maxNbColumn - 1) {
if(statesInSlice > 0) {
int symThreads = get_thread_count();
if(symThreads < 1) symThreads = 1 ;
if(symThreads > statesInSlice) symThreads = statesInSlice;
if(statesInSlice < 256) symThreads = 1 ;
std::vector<pthread_t> symWorkers(symThreads);
std::vector<SymmetrySumTask> symTasks(symThreads);
int per = (statesInSlice + symThreads -1) / symThreads ;
for(int thread = 0; thread < symThreads; thread++) {
int start = firstStateForMask + thread * per ;
int end = start + per ;
if(end > stateCount) end = stateCount;
SymmetrySumTask &task = symTasks[thread];
task.states = states;
task.partitions = partitions;
task.transitions = transitions;
task.mask_count = maskCount;
task.state_start = start;
task.state_end = end;
task.partial = 0;
pthread_create(&symWorkers[thread], NULL, accumulateSymmetryWorker, &task);
}
for(int thread = 0; thread < symThreads; thread++) {
pthread_join(symWorkers[thread], NULL);
sumSymmetric += symTasks[thread].partial ;
}
}
stateCount = firstStateForMask ;
}
}
firstStateInPrevLayer = firstStateInCurrentLayer ;
firstStateInCurrentLayer = stateCount ;
hashIndex = combinationIndexDestroy(hashIndex);
int exp2 = columns * width ;
if(columns == maxNbColumn-1) {
exp2 = (columns + 1) * width ;
result = (double)sumSymmetric / (double) (1LL << (exp2/2 )) / (double)(1LL << ((exp2+1)/2)) ;
printf("%.8f\n", result);
}
}
free(partitionBaseIndex);
free(states);
return 0 ;
}
typedef struct MaskTransition {
uint8_t surface_increment ;
uint16_t contact_mask ;
} MaskTransition ;
typedef struct MaskProfile {
int mask_union ;
int set_count ;
uint16_t part_masks[MAX_CONNECTED] ;
} MaskProfile ;
void buildPrecompute(PrecomputeData *solver_data, int width) {
solver_data->mask_count = 1 << width ;
int *bit_count_by_mask = (int*) malloc(solver_data->mask_count * sizeof(bit_count_by_mask[0]));
for(int mask = 0; mask < solver_data->mask_count; ++mask) {
int bit_count = 0 ;
for(int bit = 0; bit < width; ++bit) {
if(mask & (1 << bit)) ++bit_count ;
}
bit_count_by_mask[mask] = bit_count ;
}
solver_data->symmetry = (int*) malloc(solver_data->mask_count * sizeof(solver_data->symmetry[0]));
for(int mask = 0; mask < solver_data->mask_count; ++mask) {
int mirror = 0 ;
for(int bit = 0; bit < width; ++bit) {
mirror |= ((mask >> bit) & 1) << (width - bit - 1) ;
}
if(mirror == mask) solver_data->symmetry[mask] = 0 ;
else if(mask > mirror) solver_data->symmetry[mask] = 1;
else solver_data->symmetry[mask] = -1;
}
MaskTransition *intersections = (MaskTransition*) malloc(solver_data->mask_count * solver_data->mask_count * sizeof(intersections[0]));
for(int prev_mask = 0; prev_mask < solver_data->mask_count; ++prev_mask) {
for(int cur_mask = 0; cur_mask < solver_data->mask_count; ++cur_mask) {
int km = prev_mask & cur_mask ;
if(km==0) {
intersections[prev_mask*solver_data->mask_count + cur_mask].contact_mask = 0 ;
intersections[prev_mask*solver_data->mask_count + cur_mask].surface_increment = 0 ;
continue ;
}
for(int m = 0; m < width; m++) {
if((m>0) && (km & (1 << (m - 1))) && (cur_mask & (1 << m)) ) km |= 1 << m ;
if((m<width-1) && (km & (1 << (m + 1))) && (cur_mask & (1 << m)) ) km |= 1 << m ;
}
for(int m = width - 1; m >= 0; m--) {
if((m>0) && (km & (1 << (m - 1))) && (cur_mask & (1 << m)) ) km |= 1 << m ;
if((m<width-1) && (km & (1 << (m + 1))) && (cur_mask & (1 << m)) ) km |= 1 << m ;
}
intersections[prev_mask*solver_data->mask_count + cur_mask].surface_increment = (uint8_t) bit_count_by_mask[km] ;
intersections[prev_mask*solver_data->mask_count + cur_mask].contact_mask = (uint16_t) km ;
}
}
MaskProfile *profiles = (MaskProfile*) malloc(solver_data->mask_count * sizeof(profiles[0])) ;
int partition_slots = 0 ;
for(int mask = 0; mask < solver_data->mask_count; ++mask) {
uint16_t component_masks[MAX_CONNECTED] ;
memset(component_masks,0,sizeof(component_masks));
int component_count = 0 ;
int run = mask & 1 ;
for(int bit = 1; bit <= width; ++bit) {
if(mask & (1 << bit)) {
++run ;
} else if(run) {
component_masks[component_count++] = (uint16_t)((1 << bit) - (1 << (bit - run))) ;
run = 0 ;
}
}
assert(component_count <= width);
int partitions_by_components[MAX_CONNECTED+1] = { 1,1,2,5,14,42} ;
partition_slots += partitions_by_components[component_count] ;
profiles[mask].mask_union = mask ;
profiles[mask].set_count = component_count ;
memcpy(profiles[mask].part_masks, component_masks, sizeof(component_masks)) ;
}
solver_data->partitions = (PartitionLayout*) malloc(partition_slots * sizeof(solver_data->partitions[0]));
solver_data->first_partition_by_mask = (int*) malloc((solver_data->mask_count + 1) * sizeof(solver_data->first_partition_by_mask[0])) ;
solver_data->partition_count = 0 ;
solver_data->first_partition_by_mask[0] = 0 ;
for(int mask = 0; mask < solver_data->mask_count; ++mask) {
uint16_t component_masks[MAX_CONNECTED] ;
memcpy(component_masks, profiles[mask].part_masks, sizeof(component_masks));
int component_count = profiles[mask].set_count ;
int first = solver_data->partition_count ;
PartitionLayout *slot = nullptr ;
#define INIT_PARTITION(count) do { \
slot = solver_data->partitions + solver_data->partition_count; \
slot->id = solver_data->partition_count; \
slot->source_mask = mask; \
++solver_data->partition_count; \
memset(slot->part_masks,0,sizeof(slot->part_masks)); \
slot->set_count = (count); \
} while(0)
#define DST(i) slot->part_masks[i]
#define SRC(i) component_masks[i]
#define PART_A(i0) INIT_PARTITION(1); DST(0)=SRC(i0)
#define PART_AB(i0,i1) INIT_PARTITION(1); DST(0)=SRC(i0)|SRC(i1)
#define PART_A_B(i0,i1) INIT_PARTITION(2); DST(0)=SRC(i0); DST(1)=SRC(i1)
#define PART_ABC(i0,i1,i2) INIT_PARTITION(1); DST(0)=SRC(i0)|SRC(i1)|SRC(i2)
#define PART_AB_C(i0,i1,i2) INIT_PARTITION(2); DST(0)=SRC(i0)|SRC(i1); DST(1)=SRC(i2)
#define PART_A_B_C(i0,i1,i2) INIT_PARTITION(3); DST(0)=SRC(i0); DST(1)=SRC(i1); DST(2)=SRC(i2)
#define PART_ABCD(i0,i1,i2,i3) INIT_PARTITION(1); DST(0)=SRC(i0)|SRC(i1)|SRC(i2)|SRC(i3)
#define PART_AB_CD(i0,i1,i2,i3) INIT_PARTITION(2); DST(0)=SRC(i0)|SRC(i1); DST(1)=SRC(i2)|SRC(i3)
#define PART_ABC_D(i0,i1,i2,i3) INIT_PARTITION(2); DST(0)=SRC(i0)|SRC(i1)|SRC(i2); DST(1)=SRC(i3)
#define PART_AB_C_D(i0,i1,i2,i3) INIT_PARTITION(3); DST(0)=SRC(i0)|SRC(i1); DST(1)=SRC(i2); DST(2)=SRC(i3)
#define PART_A_B_C_D(i0,i1,i2,i3) INIT_PARTITION(4); DST(0)=SRC(i0); DST(1)=SRC(i1); DST(2)=SRC(i2); DST(3)=SRC(i3)
#define PART_ABCDE(i0,i1,i2,i3,i4) INIT_PARTITION(1); DST(0)=SRC(i0)|SRC(i1)|SRC(i2)|SRC(i3)|SRC(i4)
#define PART_ABCD_E(i0,i1,i2,i3,i4) INIT_PARTITION(2); DST(0)=SRC(i0)|SRC(i1)|SRC(i2)|SRC(i3); DST(1)=SRC(i4)
#define PART_ABC_DE(i0,i1,i2,i3,i4) INIT_PARTITION(2); DST(0)=SRC(i0)|SRC(i1)|SRC(i2); DST(1)=SRC(i3)|SRC(i4)
#define PART_ABC_D_E(i0,i1,i2,i3,i4) INIT_PARTITION(3); DST(0)=SRC(i0)|SRC(i1)|SRC(i2); DST(1)=SRC(i3); DST(2)=SRC(i4)
#define PART_AB_CD_E(i0,i1,i2,i3,i4) INIT_PARTITION(3); DST(0)=SRC(i0)|SRC(i1); DST(1)=SRC(i2)|SRC(i3); DST(2)=SRC(i4)
#define PART_AB_C_D_E(i0,i1,i2,i3,i4) INIT_PARTITION(4); DST(0)=SRC(i0)|SRC(i1); DST(1)=SRC(i2); DST(2)=SRC(i3); DST(3)=SRC(i4)
#define PART_A_B_C_D_E(i0,i1,i2,i3,i4) INIT_PARTITION(5); DST(0)=SRC(i0); DST(1)=SRC(i1); DST(2)=SRC(i2); DST(3)=SRC(i3); DST(4)=SRC(i4)
switch(component_count) {
case 0:
INIT_PARTITION(0);
break;
case 1:
PART_A(0);
break;
case 2:
PART_AB(0,1);
PART_A_B(0,1);
break;
case 3:
PART_ABC(0,1,2);
PART_AB_C(0,1,2); PART_AB_C(0,2,1); PART_AB_C(1,2,0);
PART_A_B_C(0,1,2);
break;
case 4:
PART_ABCD(0,1,2,3);
PART_ABC_D(0,1,2,3); PART_ABC_D(0,1,3,2); PART_ABC_D(0,2,3,1); PART_ABC_D(1,2,3,0);
PART_AB_CD(0,1,2,3); PART_AB_CD(0,3,1,2);
PART_AB_C_D(0,1,2,3); PART_AB_C_D(0,2,1,3); PART_AB_C_D(0,3,1,2); PART_AB_C_D(1,2,0,3); PART_AB_C_D(1,3,0,2); PART_AB_C_D(2,3,0,1);
PART_A_B_C_D(0,1,2,3);
break;
case 5:
PART_ABCDE(0,1,2,3,4);
PART_ABCD_E(0,1,2,3,4); PART_ABCD_E(0,1,2,4,3); PART_ABCD_E(0,1,3,4,2); PART_ABCD_E(0,2,3,4,1); PART_ABCD_E(1,2,3,4,0);
PART_ABC_DE(0,1,2,3,4); PART_ABC_DE(0,1,4,2,3); PART_ABC_DE(0,3,4,1,2); PART_ABC_DE(2,3,4,0,1); PART_ABC_DE(1,2,3,0,4);
PART_ABC_D_E(0,1,2,3,4); PART_ABC_D_E(0,1,3,2,4); PART_ABC_D_E(0,2,3,1,4); PART_ABC_D_E(1,2,3,0,4); PART_ABC_D_E(0,1,4,2,3);
PART_ABC_D_E(0,2,4,1,3); PART_ABC_D_E(1,2,4,0,3); PART_ABC_D_E(0,3,4,1,2); PART_ABC_D_E(1,3,4,0,2); PART_ABC_D_E(2,3,4,0,1);
PART_AB_CD_E(0,1,2,3,4); PART_AB_CD_E(0,3,1,2,4); PART_AB_CD_E(0,1,2,4,3); PART_AB_CD_E(0,4,1,2,3); PART_AB_CD_E(0,1,3,4,2);
PART_AB_CD_E(0,4,1,3,2); PART_AB_CD_E(0,2,3,4,1); PART_AB_CD_E(0,4,2,3,1); PART_AB_CD_E(1,2,3,4,0); PART_AB_CD_E(1,4,2,3,0);
PART_AB_C_D_E(0,1,2,3,4); PART_AB_C_D_E(0,2,1,3,4); PART_AB_C_D_E(0,3,1,2,4); PART_AB_C_D_E(0,4,1,2,3); PART_AB_C_D_E(1,2,0,3,4);
PART_AB_C_D_E(1,3,0,2,4); PART_AB_C_D_E(1,4,0,2,3); PART_AB_C_D_E(2,3,0,1,4); PART_AB_C_D_E(2,4,0,1,3); PART_AB_C_D_E(3,4,0,1,2);
PART_A_B_C_D_E(0,1,2,3,4);
break;
}
solver_data->first_partition_by_mask[mask+1] = solver_data->partition_count ;
for(int ip = first; ip < solver_data->partition_count; ++ip) {
int count = solver_data->partitions[ip].set_count;
for(int i = 1; i < count; ++i) {
for(int j = i; j < count; ++j) {
if(solver_data->partitions[ip].part_masks[i-1] > solver_data->partitions[ip].part_masks[j]) {
uint16_t tmp = solver_data->partitions[ip].part_masks[i-1] ;
solver_data->partitions[ip].part_masks[i-1] = solver_data->partitions[ip].part_masks[j] ;
solver_data->partitions[ip].part_masks[j] = tmp ;
}
}
}
}
}
solver_data->transitions = (PartitionTransition*) calloc(solver_data->partition_count * solver_data->mask_count, sizeof(solver_data->transitions[0]));
for(int partition_id = 0; partition_id < solver_data->partition_count; ++partition_id) {
const PartitionLayout *from_partition = solver_data->partitions + partition_id ;
for(int mask = 0; mask < solver_data->mask_count; ++mask) {
int nbSet = 0 ;
uint16_t merged[MAX_CONNECTED] ;
memset(merged,0,sizeof(merged));
int global_contact = 0 ;
for(int i = 0; i < from_partition->set_count; ++i) {
uint16_t new_contact = intersections[from_partition->part_masks[i] * solver_data->mask_count + mask].contact_mask ;
if(new_contact) {
for(int j = 0; j < nbSet; ++j) {
if((new_contact & merged[j])) {
new_contact = intersections[(new_contact | merged[j]) * solver_data->mask_count + mask].contact_mask ;
for(int k = j + 1; k < nbSet; ++k) merged[k - 1] = merged[k] ;
--nbSet ;
--j ;
merged[nbSet] = 0 ;
}
}
int pos = nbSet - 1 ;
while(pos >= 0 && merged[pos] > new_contact) {
merged[pos + 1] = merged[pos] ;
--pos ;
}
merged[pos + 1] = new_contact ;
++nbSet ;
}
global_contact |= new_contact ;
}
if(global_contact != mask) {
global_contact ^= mask ;
for(int k = 0; k < profiles[mask].set_count; ++k) {
if((profiles[mask].part_masks[k] & global_contact) == profiles[mask].part_masks[k] ) {
int pos = nbSet - 1 ;
while(pos >= 0 && merged[pos] > profiles[mask].part_masks[k]) {
merged[pos + 1] = merged[pos] ;
--pos ;
}
merged[pos + 1] = profiles[mask].part_masks[k] ;
++nbSet ;
}
}
}
int target_partition = solver_data->first_partition_by_mask[mask] ;
while(target_partition < solver_data->first_partition_by_mask[mask + 1] && memcmp(merged, solver_data->partitions[target_partition].part_masks, sizeof(merged)) != 0) {
++target_partition ;
}
assert(target_partition != solver_data->first_partition_by_mask[mask + 1]) ;
PartitionTransition *transition = solver_data->transitions + partition_id * solver_data->mask_count + mask ;
transition->next_partition = (int)target_partition ;
for(int id = 0; id < solver_data->partitions[target_partition].set_count; ++id) {
transition->added_surface[id] = (uint8_t) bit_count_by_mask[solver_data->partitions[target_partition].part_masks[id]] ;
}
for(int i = 0; i < from_partition->set_count; ++i) {
if(from_partition->part_masks[i] & mask) {
for(int id = 0; id < solver_data->partitions[target_partition].set_count; ++id){
if(from_partition->part_masks[i] & solver_data->partitions[target_partition].part_masks[id]) {
transition->source_index[i] = (uint8_t)(id + 1) ;
break ;
}
assert(id < solver_data->partitions[target_partition].set_count) ;
}
}
}
}
}
free(profiles);
free(bit_count_by_mask);
free(intersections);
}
CombinationIndex * combinationIndexAlloc(IndexValue N, IndexValue k) {
if(k<2) return NULL ;
CombinationIndex * NS = (CombinationIndex*) calloc(1, sizeof(NS[0])) ;
NS->N = N+1 ;
NS->k = k ;
NS->index = (IndexValue*) malloc((k-1) * (N+1) * sizeof(NS->index[0]));
IndexValue *ids = NS->index+ ((k-2)*(N+1) );
ids[0] = 0 ;
for(IndexValue in=1;in<=N;in++){
ids[in] = ids[in-1] + N + 2 - in ;
}
for(int ik=k-3;ik>=0;ik--) {
IndexValue *idsAnt = ids ;
ids = NS->index+ (ik*(N+1));
ids[0] = 0 ;
IndexValue indAnt = idsAnt[N]+1 ;
for(IndexValue in=1;in<=N;in++){
ids[in] = ids[in-1] + indAnt - idsAnt[in-1] ;
}
}
return NS ;
}
CombinationIndex * combinationIndexDestroy(CombinationIndex * NS) {
if(NS != NULL) {
free(NS->index);
free(NS);
}
return NULL ;
}
IndexValue combinationIndexSize(const CombinationIndex *NS,int ks) {
if(ks<=1) {
return ks ? NS->N : 1 ;
}
return NS->index[(NS->k-ks)*(NS->N) + (NS->N - 1)] + 1 ;
}
IndexValue combinationIndexRank(const CombinationIndex *NS,int ks,SurfaceValue *sum) {
int k = NS->k ;
int N = NS->N ;
IndexValue * index = NS->index + (k-ks)*N ;
switch(ks) {
case 0: return 0;
case 1: return sum[0] ;
case 2: return index[sum[0]]+sum[1] ;
case 3: return index[sum[0]] + index[N+sum[0]+sum[1]]-index[N+sum[0]]+sum[2];
case 4: return index[sum[0]] + index[N+sum[0]+sum[1]]-index[N+sum[0]]
+ index[2*N+sum[0]+sum[1]+sum[2]]-index[2*N+sum[0]+sum[1]] +sum[3];
default : {
IndexValue ksum = *sum++ ;
IndexValue ind = index[ksum] ;
for(int ik=1;ik< ks-1;ik++) {
IndexValue ksumNext = ksum + *sum++ ;
ind += index[N*ik+ksumNext]-index[N*ik+ksum];
ksum=ksumNext ;
}
return ind+ *sum ;
}
}
}
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.nio.file.*;
import java.util.*;
import java.util.regex.*;
public class Euler701 {
private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");
private static String parseOutput(String stdout) {
String[] lines = stdout.split("\\R");
List<String> nonEmpty = new ArrayList<>();
for (String line : lines) {
String t = line.trim();
if (!t.isEmpty()) {
nonEmpty.add(t);
}
}
if (nonEmpty.isEmpty()) {
return "";
}
List<String> answers = new ArrayList<>();
List<String> equals = new ArrayList<>();
for (String line : nonEmpty) {
Matcher m1 = ANSWER_RE.matcher(line);
if (m1.find()) {
answers.add(m1.group(1).trim());
}
Matcher m2 = EQUAL_RE.matcher(line);
if (m2.find()) {
equals.add(m2.group(1).trim());
}
}
if (!answers.isEmpty()) {
return answers.get(answers.size() - 1);
}
if (!equals.isEmpty()) {
return equals.get(equals.size() - 1);
}
return nonEmpty.get(nonEmpty.size() - 1);
}
private static String pickCompiler() throws Exception {
for (String compiler : List.of("clang++", "g++")) {
Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
.redirectErrorStream(true)
.start();
String out = new String(probe.getInputStream().readAllBytes());
int rc = probe.waitFor();
if (rc == 0 && !out.trim().isEmpty()) {
return compiler;
}
}
throw new RuntimeException("No C++ compiler found (clang++/g++).");
}
private static Path cppSource(Path root) {
return root.resolve("solutionsCpp").resolve("Euler701.cpp");
}
private static boolean shouldSkipCheckpoints(Path root) {
Path src = cppSource(root);
try {
String text = Files.readString(src);
return text.contains("--skip-checkpoints");
} catch (Exception ex) {
return false;
}
}
private static Path ensureBridgeBinary() throws Exception {
Path root = Paths.get(System.getProperty("user.dir"));
Path src = cppSource(root);
Path bin = root.resolve("solutionsCpp").resolve(".euler701_java_bridge");
boolean rebuild = Files.notExists(bin)
|| Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;
if (rebuild) {
String compiler = pickCompiler();
Process compile = new ProcessBuilder(
compiler,
"-std=c++17",
"-O2",
src.toString(),
"-o",
bin.toString())
.inheritIO()
.start();
if (compile.waitFor() != 0) {
throw new RuntimeException("Failed to compile Euler701 C++ bridge.");
}
}
return bin;
}
private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
List<String> cmd = new ArrayList<>();
cmd.add(bin.toString());
if (shouldSkipCheckpoints(root)) {
cmd.add("--skip-checkpoints");
}
Process first = new ProcessBuilder(cmd)
.directory(root.toFile())
.redirectErrorStream(true)
.start();
String out = new String(first.getInputStream().readAllBytes());
int rc = first.waitFor();
if (rc == 0) {
return out;
}
Process second = new ProcessBuilder(cmd)
.directory(srcDir.toFile())
.redirectErrorStream(true)
.start();
String out2 = new String(second.getInputStream().readAllBytes());
int rc2 = second.waitFor();
if (rc2 == 0) {
return out2;
}
throw new RuntimeException("Euler701 C++ bridge failed.\n" + out + "\n" + out2);
}
private static String solveViaCppBridge() throws Exception {
Path root = Paths.get(System.getProperty("user.dir"));
Path src = cppSource(root);
Path bin = ensureBridgeBinary();
String out = runBridge(bin, root, src.getParent());
String parsed = parseOutput(out);
if (parsed.isEmpty()) {
throw new RuntimeException("Euler701 C++ bridge produced empty output.");
}
return parsed;
}
public static void main(String[] args) throws Exception {
System.out.println(solveViaCppBridge());
}
}