Problem 300: Protein Folding

Project Euler Problem 300 Solution

EulerSolve provides an optimized solution for Project Euler Problem 300, Protein Folding, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every binary HP sequence of length $n$, we place the chain on the square lattice as a self-avoiding walk and ask for the maximum possible number of H-H contacts. Project Euler asks for the average of that maximum over all $2^n$ sequences. The code solves the case $$n=15.$$ Mathematical Approach 1. Folding shapes are self-avoiding walks A fold is a self-avoiding walk $$p_0,p_1,\dots,p_{n-1}\in\mathbb Z^2$$ with unit steps and no repeated lattice site. Because translation and rotation do not change contact counts, the implementation fixes $$p_0=(0,0),\qquad p_1=(1,0),$$ and then performs a DFS over the remaining $n-2$ steps. This removes the large symmetry factor before enumeration even starts. 2. A fold matters only through its contact map For a finished walk, the only geometric information relevant to scoring is which non-consecutive monomer indices end up adjacent on the lattice. For $i<j$ with $j\ge i+2$, define $$C_{ij}=1 \iff |p_i-p_j|_1=1.$$ This produces a contact map, which the code stores as a bitset over all pairs $(i,j)$ with $j\ge i+2$. Different lattice walks can produce the same contact map, and then they have exactly the same score for every HP sequence. So the geometry phase ends by deduplicating maps with hashing....

Detailed mathematical approach

Problem Summary

A protein is modelled as a string $s=s_0s_1\cdots s_{n-1}$ of elements $\mathrm H$ (hydrophobic) and $\mathrm P$ (polar). Orientation matters: HPP and PPH are different strings, so there are $2^n$ strings of length $n$. A folding puts the elements on points of the square lattice $\mathbb Z^2$ so that consecutive elements sit on neighbouring points and no point is used twice. An H-H contact point is a pair of $\mathrm H$ elements on neighbouring points, consecutive ones included, and an optimal folding of a string has as many of them as possible. With $\mathrm H$ and $\mathrm P$ equally likely at every position, all $2^n$ strings are equally likely, and we must find the exact average number of H-H contact points in an optimal folding of a random string of length $n=15$. The statement shows HHPPHHHPHHPH in a folding with $6$ contact points and in an optimal one with $9$, and it gives the average $850/2^8=3.3203125$ for $n=8$. ^[1]

The whole solution in one sentence: list the foldings of $15$ elements up to rotation and translation, keep the distinct sets of non-consecutive neighbour pairs that they produce, score each set against all $2^{15}$ strings at once with a bitmask recurrence, and average the best scores after adding the consecutive H-H pairs, which every folding has.

Mathematical Approach

Fifteen elements have few enough foldings to list them all; the work lies in sharing that list between all $32768$ strings.

1. Foldings and what is counted

A folding of length $n$ is a self-avoiding walk $w=(p_0,\dots,p_{n-1})$ on $\mathbb Z^2$: consecutive points are at distance $1$, and no point is visited twice. ^[3] Two lattice points are neighbours when one coordinate differs by $1$ and the other is equal. For a string $s$, let $E(s,w)$ be the number of pairs $i\lt j$ with $s_i=s_j=\mathrm H$ whose points are neighbours. The pairs with $j=i+1$, the bonds, are neighbours in every folding, and the statement counts them: HHPPHHHPHHPH has $4$ bonds between two $\mathrm H$, and listing all its foldings gives totals from $4$ to $9$, so its optimum $9$ consists of the $4$ bonds and $5$ further contacts. The standard HP model of Dill rewards only the H-H pairs that are not bonded. ^[2] The difference is a number fixed by the string, so both counts have the same optimal foldings, and the checkpoint confirms the count with bonds: the $256$ strings of length $8$ give a total of $850$ with bonds and $402$ without. We write

$$\begin{aligned}F(s)&=\max_{w}E(s,w),\\ A(n)&=\frac{1}{2^n}\sum_{s}F(s)\end{aligned}$$

for the optimum of a string and the average; we must compute $A(15)$.

2. The bonds contribute a fixed amount

Let $b(s)$ be the number of bonds whose two elements are both $\mathrm H$, and let $G(s)$ be the largest number of non-bonded H-H contacts in a folding of $s$. The bonds are present in every folding, so

$$F(s)=b(s)+G(s).$$

Each of the $n-1$ bonds joins two $\mathrm H$ in exactly $2^{n-2}$ strings, because its two letters are fixed and the other $n-2$ are free. Hence

$$\sum_{s}b(s)=(n-1)\,2^{n-2},$$

and the bonds add exactly $(n-1)/4$ to the average, for example $448/256=1.75$ for $n=8$. All the difficulty lies in $G$.

3. Only the contact map matters

The contact map $C(w)$ of a folding is the set of pairs $(i,j)$ with $j\ge i+2$ whose points are neighbours; for example, the square folding $(0,0)$, $(1,0)$, $(1,1)$, $(0,1)$ of four elements has $C(w)=\{(0,3)\}$. For any set $C$ of such pairs, let $\sigma_C(s)$ be the number of pairs in $C$ whose two elements are both $\mathrm H$ in $s$. The non-bonded H-H contacts of $s$ in $w$ are exactly the pairs counted by $\sigma_{C(w)}(s)$, so they depend on $w$ only through its contact map, and

$$G(s)=\max_{C}\,\sigma_C(s),$$

where $C$ runs over the distinct contact maps of all foldings of length $n$. This is the decisive compression: the $593\,611$ walks searched for $n=15$ produce only $12495$ distinct maps, and each map is compared with the strings only once.

4. A chessboard colouring

Colour every lattice point $(x,y)$ by the parity of $x+y$. A unit step changes $x+y$ by $\pm1$, so neighbours always have different colours: the square grid is a bipartite graph. ^[6] Along a walk the colours alternate, so $p_i$ has the colour of $p_0$ exactly when $i$ is even. Two elements $i\lt j$ can therefore be neighbours only when $j-i$ is odd, and a non-bonded contact needs $j-i\ge3$. For $n=15$ the candidates are the pairs with $j-i\in\{3,5,\dots,13\}$, $12+10+8+6+4+2=42$ of them, and each of them occurs in some folding. Counting neighbours also bounds the size of a map: a point has four neighbours, an inner element needs two of them for its bonds and an end element needs one, so counting every contact from both of its elements gives

$$|C(w)|\le\frac{3+3+2(n-2)}{2}=n+1.$$

For $n=15$ the largest map that occurs has $8$ pairs.

5. Listing the foldings

Translations, rotations by multiples of $90^\circ$ and reflections map neighbours to neighbours, so they do not change contact maps. A translation moves $p_0$ to $(0,0)$, and then exactly one rotation turns the first step into the step to $(1,0)$. Every folding is therefore carried by a translation and a rotation onto exactly one walk with $p_0=(0,0)$ and $p_1=(1,0)$, and this walk has the same contact map. If $c_k$ is the number of $k$-step self-avoiding walks from the origin, there are $c_{n-1}/4$ such walks: $c_7=2172$ gives $543$ walks for $n=8$, and $c_{14}=2\,374\,444$ gives $593\,611$ for $n=15$. ^[4] These numbers grow roughly like $\mu^k$ with the connective constant $\mu\approx2.638$ of the square lattice, ^[5] so a complete list is possible only for short chains. The search does not exploit the reflection $(x,y)\mapsto(x,-y)$, which pairs every walk except the straight one with a different walk of the same map, and these duplicates are removed like all others.

The walks are produced by a depth-first search with backtracking: place $p_{k+1}$ on each free neighbour of $p_k$ in turn, continue from there, and free the point again afterwards. A map is recorded only when all $n$ points are placed. Since adding points can only add contacts, it is tempting to score partial walks as well, but a partial walk can be a dead end whose last point has no free neighbour, and then it has no completion at all. ^[3] The first dead ends appear at $8$ points, and the search for $n=15$ meets $2228$ of them. A direct search that counted the H-H pairs among the placed elements, bonds included, and maximised over partial walks as well would report the total $1999$ instead of the correct $1997$ for length $9$.

6. Worked example: $n=4$, $5$ and $6$

For $n=4$ the $9$ walks give two maps: $7$ walks have no contact, and the $2$ walks around a unit square give $\{(0,3)\}$. So $G(s)=1$ exactly when $s_0=s_3=\mathrm H$, which holds for $4$ strings, and with the $3\cdot2^2=12$ bonds the total is $16$ and $A(4)=1$.

For $n=5$ the $25$ walks give three maps: the empty map from $17$ walks, $\{(0,3)\}$ from $4$ and $\{(1,4)\}$ from $4$. The two contacts never occur together: $(0,3)$ makes $p_0p_1p_2p_3$ a unit square, so $p_1$ and $p_3$ are diagonal, their only common neighbours $p_0$ and $p_2$ are taken, and $p_4$ cannot touch both of them. Hence $G(s)=1$ when $s_0=s_3=\mathrm H$ or $s_1=s_4=\mathrm H$, and $G(s)=0$ otherwise. Each condition holds for $8$ strings and both hold for $2$, so $G$ sums to $8+8-2=14$. The bonds give $4\cdot2^3=32$, and

$$A(5)=\frac{32+14}{32}=\frac{23}{16}=1.4375.$$

For $n=6$ the $71$ walks give $8$ maps: the empty map, the single contacts $(0,3)$, $(1,4)$ and $(2,5)$, and the pairs $\{(0,3),(2,5)\}$, $\{(0,3),(0,5)\}$, $\{(0,5),(1,4)\}$ and $\{(0,5),(2,5)\}$. The contact $(0,5)$ never comes alone: it closes the chain into a cycle of six points, which must bound a $1\times2$ rectangle, and the midpoints of the long sides form a second contact. A complete enumeration gives the totals $\sum_sF(s)$ equal to $1$, $4$, $16$, $46$, $133$, $335$ and $850$ for $n=2,\dots,8$.

7. Scoring one map against all strings

Encode a string by the set $m\subseteq\{0,\dots,n-1\}$ of its $\mathrm H$ positions, stored as a bitmask whose bit $i$ is set when $s_i=\mathrm H$. For a map $C$, let $N_C(i)$ be the set of contact partners of $i$, and write $\sigma_C(m)$ for the number of pairs of $C$ inside $m$, so that $\sigma_C(\varnothing)=0$. If $i$ is the smallest element of a nonempty set $m$ and $m'=m\setminus\{i\}$, then

$$\sigma_C(m)=\sigma_C(m')+\bigl|N_C(i)\cap m'\bigr|.$$

Indeed, a pair of $C$ inside $m$ either avoids $i$, and then it lies inside $m'$, or it consists of $i$ and one of its partners in $m'$. As a number, $m'=m-2^i$ is smaller than $m$, so filling a table in increasing order of $m$ makes $\sigma_C(m')$ available when it is needed, and every entry costs one step. This is dynamic programming over subsets, and the smallest element is the lowest set bit, which the expression m & -m isolates. ^[7] For example, take the folding of $n=6$ around the $1\times2$ rectangle, with map $\{(0,5),(1,4)\}$, and the string HHPPHH, so $m=\{0,1,4,5\}$. Removing $0$, $1$, $4$ and $5$ in turn gives, from the innermost value outwards,

$$\begin{aligned}\sigma(\{5\})&=0,\\ \sigma(\{4,5\})&=0,\\ \sigma(\{1,4,5\})&=1,\\ \sigma(\{0,1,4,5\})&=2.\end{aligned}$$

With its $2$ bonds the string has $4$ contacts in this folding, and the enumeration confirms that $4$ is its optimum.

8. The best map for every string

Let $B(m)$ be the maximum of $\sigma_C(m)$ over the distinct maps. The search of Section 5 meets every folding up to translation and rotation, and every map that it records belongs to a real folding, so $B(m)=G(s)$ for the string $s$ whose $\mathrm H$ set is $m$. The bonds of $s$ are the positions $i$ with $i\in m$ and $i+1\in m$; bit $i$ of m & (m >> 1) is set exactly when bits $i$ and $i+1$ of $m$ are, so $b(s)$, which we also write $b(m)$, is the number of set bits of this mask. Together,

$$\begin{aligned}A(n)&=\frac{T_n}{2^n},\\ T_n&=\sum_{m=0}^{2^n-1}\bigl(B(m)+b(m)\bigr).\end{aligned}$$

For $n=8$ the search visits $543$ walks and finds $41$ maps, the best scores add up to $402$ and the bonds to $7\cdot2^6=448$, and $T_8=850$ gives the value $A(8)=3.3203125$ of the statement.

9. An exact decimal

The denominator is a power of two, so the average has a terminating decimal expansion, like every fraction whose denominator has no prime factors other than $2$ and $5$. ^[9] More precisely, write $T_n/2^n=t/2^k$ in lowest terms. If $k\ge1$, then $t$ is odd and

$$\frac{t}{2^k}=\frac{5^k\,t}{10^k},$$

where $5^kt$ is odd, so its last digit is not $0$ and the expansion has exactly $k\le n$ digits after the decimal point. Long division produces them one at a time: after each digit the remainder $r$ becomes $10r\bmod 2^n$, so after $d$ digits it equals $10^d r_0\bmod 2^n$, which vanishes once $d\ge n$. For $n=8$, $850/256=425/2^7$ has the $7$ decimals of $3.3203125$.

10. Shortcuts the code does not need

If $C\subseteq C'$, then $\sigma_C(m)\le\sigma_{C'}(m)$ for every $m$, so only the maps that are maximal under inclusion matter, although the code keeps all of them: $20$ of the $41$ maps for $n=8$, and $3409$ of the $12495$ for $n=15$. Reversing the chain, $i\mapsto n-1-i$, shows that a string and its reverse have the same optimum; both still count in the average. In general, finding an optimal folding in the HP model is NP-hard even on the two-dimensional square lattice, ^[2] so no fast exact method is expected for long chains; for $n=15$, exhaustive search is cheap.

How the Code Works

The C++ program ^[10] starts in main: parse_arguments reads --skip-checkpoints and --n= followed by digits, requires the value to lie from $2$ to $20$ (default $15$) and makes main exit with code $1$ when parsing fails; unless skipped, run_checkpoints requires the string 3.3203125 for $n=8$, with exit code $2$ otherwise; then solve_average_decimal builds a FoldingSolver and returns the exact decimal, which main prints.

The constructor of FoldingSolver allocates the flags occupied for a grid of side $2n+1$ whose centre $(n,n)$ will hold $p_0$; every point of a walk lies at most $n-1$ steps from $p_0$, so the walk and the neighbours that it tests stay inside the grid. It lists in pair_list all pairs $(i,j)$ with $j\ge i+2$ in lexicographic order, $91$ pairs for $n=15$, and the position $p$ of a pair in this list becomes its bit number. All are tested, including the $49$ that Section 4 rules out, whose bits stay $0$.

enumerate_maps marks $p_0$ and $p_1=p_0+(1,0)$ as occupied and calls dfs(1), the search of Section 5. A call dfs(idx) tries the directions $+x$, $-x$, $+y$ and $-y$ from the last placed point, skips occupied points, marks the new point, recurses and clears the mark again. When idx reaches $n-1$, add_current_contact_map tests every pair of pair_list for Manhattan distance $1$ and sets bit $p$ of a ContactKey, two 64-bit words lo and hi for the bits $0$ to $63$ and $64$ to $127$. The keys go into an std::unordered_set, which discards repeated maps (Section 3); its hash ContactKeyHash multiplies the two words by large odd constants and combines the products.

solve_numerator converts every key into an array adj of 16-bit masks, where bit $j$ of adj[i] is set when $i$ and $j$ form a contact: the sets $N_C(i)$ of Section 7. Two byte arrays of length $2^n$ hold score, the table $\sigma_C$ of the current map, and best, the running maximum $B$ of Section 8. For each map, the loop over mask from $1$ to $2^n-1$ computes lb = mask & -mask, the index i = __builtin_ctz(lb) of the lowest set bit, prev = mask ^ lb and the increment __builtin_popcount(prev & adj[i]), exactly the recurrence of Section 7. ^[8] __builtin_ctz is undefined for $0$, which cannot occur because mask is at least $1$, and a byte holds every score because a map has at most $n+1=16$ pairs. Finally a 64-bit sum of best[mask] and __builtin_popcount(mask & (mask >> 1)) over all masks gives $T_n$, and exact_average_decimal performs the long division of Section 9 in 128-bit integers. These widths cover $n=15$, which needs $15$ mask bits and $91$ key bits; the 16-bit masks suffice only for $n\le16$, and the 128 key bits hold the $(n-1)(n-2)/2$ pairs only for $n\le17$, although the argument check accepts values up to $20$.

The Python version is a single function solve with $n=15$ fixed and no checkpoint. Its keys are Python integers, which have unlimited precision, so one number holds a whole map, and a set collects them; bytearray objects serve as occupied, as best and as a fresh score table for each map, the lowest set bit has index lb.bit_length() - 1, and a table built by popcount[m] = popcount[m >> 1] + (m & 1) counts set bits. It raises the recursion limit, although the search is only $n-1$ levels deep. The Java version follows the C++ structure, also with $n=15$ fixed and no checkpoint: ContactKey implements equals and a multiplicative hashCode, the masks are int arrays, the tables are byte arrays, Integer.numberOfTrailingZeros and Integer.bitCount replace the GCC built-ins, and the long division uses long.

Complexity and Verification

For $n=15$ the search makes $942\,710$ calls of dfs, one for each partial walk of $2$ to $15$ points, and reaches $593\,611$ complete walks, each of which tests $91$ pairs, about $5.4\cdot10^7$ tests in all. The scoring makes $M(2^n-1)$ recurrence steps for $M$ distinct maps, here $12495\cdot32767=409\,423\,665$ steps of a few bit operations each, and it dominates. In total the time is $O(n^2c_{n-1}+M\,2^n)$ and the memory $O(nM+2^n)$. In the C++ program the data are $12495$ keys of $16$ bytes, kept in the hash set and in a copied vector, as many mask arrays of $40$ bytes and two tables of $32768$ bytes, and the whole process stays below $4$ MB. On an Apple-silicon Mac (clang with -O2, OpenJDK 17, CPython 3.11) the C++ program takes about half a second, Java about a second and Python about forty seconds.

The C++ program checks itself with run_checkpoints for $n=8$; the Python and Java versions contain no check. For this article the following was verified separately, outside the published solutions.

An independent brute force searches the walks separately for every string, counts its H-H contacts, bonds included, as the walk grows and keeps the maximum over complete walks; it uses no contact maps and no recurrence. It agrees with the C++ program for every $n$ from $2$ to $13$, and its totals for $n\le9$ are those of Sections 5 and 6.
A second computation for $n=15$ listed only the walks whose first turn, if any, goes upward, $296\,806$ of them, found the same $12495$ maps, kept the $3409$ maximal ones and counted the H-H pairs of every string directly. Its output is identical to that of the C++ program.
The C++, Java and Python programs print the same value for $n=15$.
The walk counts $543$ and $593\,611$ agree with $c_7/4$ and $c_{14}/4$ from the OEIS. ^[4]

Footnotes and References

The references below supply the background facts used above. The split into bonds and contact maps, the subset recurrence, the small cases and all counts are derived or computed explicitly in this article.

Project Euler 300 — Protein Folding. The official statement describes strings of H and P elements whose orientation matters, folded in two dimensions; it shows HHPPHHHPHHPH folded with 6 and, optimally, with 9 H-H contact points, gives 850/2⁸ = 3.3203125 for length 8 and asks for the exact average for length 15.
Wikipedia — Hydrophobic-polar protein folding model. Proposed by Ken Dill in 1985, the model folds a chain as a self-avoiding walk on a lattice and rewards adjacent H residues that are not covalently bound; finding an optimal folding is NP-hard on the two-dimensional square lattice (Crescenzi et al., 1998).
Wikipedia — Self-avoiding walk. A self-avoiding walk is a sequence of moves on a lattice that never visits the same point twice; a growing walk is trapped when its end has no empty neighbouring site.
OEIS A001411 — Number of n-step self-avoiding walks on square lattice. The terms a(7) = 2172 and a(14) = 2374444 give the walk counts 543 and 593611 with a fixed first step.
Wikipedia — Connective constant. The connective constant μ = lim cₙ^(1/n) is about 2.63815853 for the square lattice, and cₙ is conjectured to grow like μⁿ times a power of n.
Wikipedia — Lattice graph. The square grid graph joins two integer points at distance one, and all square grid graphs are bipartite.
cp-algorithms — Bit manipulation. An integer can record a set of objects, one bit each; n & −n extracts the last set bit, and __builtin_popcount and __builtin_ctz count set bits and trailing zeros.
GCC manual — Bit Operation Builtins. __builtin_ctz returns the number of trailing 0-bits and is undefined for 0; __builtin_popcount returns the number of 1-bits.
Wikipedia — Repeating decimal. A fraction in lowest terms has a non-terminating, repeating expansion exactly when its denominator has a prime factor other than 2 and 5.
C++ — Euler300.cpp. The linked immutable C++ revision contains FoldingSolver with dfs, add_current_contact_map, enumerate_maps, solve_numerator and exact_average_decimal, and run_checkpoints with the checkpoint for $n=8$. It is the implementation record for this article.

Problem 300 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <thread>
#include <unordered_set>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

struct Options {
    int n = 15;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 2 && options.n <= 20;
}

struct ContactKey {
    u64 lo;
    u64 hi;

    bool operator==(const ContactKey& other) const {
        return lo == other.lo && hi == other.hi;
    }
};

struct ContactKeyHash {
    std::size_t operator()(const ContactKey& k) const {
        const u64 x = k.lo * 11400714819323198485ull;
        const u64 y = k.hi * 14029467366897019727ull;
        return static_cast<std::size_t>(x ^ (y >> 1));
    }
};

struct FoldingSolver {
    int n;
    int grid_side;
    int offset;

    std::vector<int> px;
    std::vector<int> py;
    std::vector<unsigned char> occupied;

    std::array<int, 4> dx{1, -1, 0, 0};
    std::array<int, 4> dy{0, 0, 1, -1};

    std::vector<std::pair<int, int>> pair_list;
    std::vector<std::vector<int>> pair_index;

    std::unordered_set<ContactKey, ContactKeyHash> contact_maps;

    explicit FoldingSolver(int n_) : n(n_) {
        grid_side = 2 * n + 1;
        offset = n;
        px.assign(static_cast<std::size_t>(n), 0);
        py.assign(static_cast<std::size_t>(n), 0);
        occupied.assign(static_cast<std::size_t>(grid_side * grid_side), 0);

        pair_index.assign(static_cast<std::size_t>(n), std::vector<int>(static_cast<std::size_t>(n), -1));
        int idx = 0;
        for (int i = 0; i < n; ++i) {
            for (int j = i + 2; j < n; ++j) {
                pair_index[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] = idx;
                pair_list.push_back({i, j});
                ++idx;
            }
        }

        contact_maps.reserve(1 << 16);
    }

    int cell_id(int x, int y) const {
        return y * grid_side + x;
    }

    void add_current_contact_map() {
        ContactKey key{0, 0};

        for (int p = 0; p < static_cast<int>(pair_list.size()); ++p) {
            const int i = pair_list[static_cast<std::size_t>(p)].first;
            const int j = pair_list[static_cast<std::size_t>(p)].second;
            const int manhattan = std::abs(px[static_cast<std::size_t>(i)] - px[static_cast<std::size_t>(j)]) +
                                  std::abs(py[static_cast<std::size_t>(i)] - py[static_cast<std::size_t>(j)]);
            if (manhattan != 1) {
                continue;
            }

            if (p < 64) {
                key.lo |= (1ULL << p);
            } else {
                key.hi |= (1ULL << (p - 64));
            }
        }

        contact_maps.insert(key);
    }

    void dfs(int idx) {
        if (idx == n - 1) {
            add_current_contact_map();
            return;
        }

        const int cx = px[static_cast<std::size_t>(idx)];
        const int cy = py[static_cast<std::size_t>(idx)];

        for (int dir = 0; dir < 4; ++dir) {
            const int nx = cx + dx[static_cast<std::size_t>(dir)];
            const int ny = cy + dy[static_cast<std::size_t>(dir)];
            const int cid = cell_id(nx, ny);
            if (occupied[static_cast<std::size_t>(cid)] != 0) {
                continue;
            }
            occupied[static_cast<std::size_t>(cid)] = 1;
            px[static_cast<std::size_t>(idx + 1)] = nx;
            py[static_cast<std::size_t>(idx + 1)] = ny;
            dfs(idx + 1);
            occupied[static_cast<std::size_t>(cid)] = 0;
        }
    }

    std::vector<ContactKey> enumerate_maps() {
        const int x0 = offset;
        const int y0 = offset;
        const int x1 = offset + 1;
        const int y1 = offset;

        occupied[static_cast<std::size_t>(cell_id(x0, y0))] = 1;
        occupied[static_cast<std::size_t>(cell_id(x1, y1))] = 1;

        px[0] = x0;
        py[0] = y0;
        px[1] = x1;
        py[1] = y1;

        dfs(1);

        std::vector<ContactKey> maps;
        maps.reserve(contact_maps.size());
        for (const ContactKey& key : contact_maps) {
            maps.push_back(key);
        }
        return maps;
    }

    static std::string exact_average_decimal(u64 numerator, int n) {
        const u128 denominator = static_cast<u128>(1) << n;
        const u128 num = numerator;
        const u128 integer_part = num / denominator;
        u128 rem = num % denominator;

        std::string out = std::to_string(static_cast<u64>(integer_part));
        if (rem == 0) {
            return out;
        }

        out.push_back('.');
        while (rem != 0) {
            rem *= 10;
            const u64 digit = static_cast<u64>(rem / denominator);
            out.push_back(static_cast<char>('0' + digit));
            rem %= denominator;
        }
        return out;
    }

    u64 solve_numerator() {
        const std::vector<ContactKey> maps = enumerate_maps();

        std::vector<std::array<std::uint16_t, 20>> adj_masks;
        adj_masks.reserve(maps.size());

        for (const ContactKey& key : maps) {
            std::array<std::uint16_t, 20> adj{};
            adj.fill(0);

            for (int p = 0; p < static_cast<int>(pair_list.size()); ++p) {
                bool on = false;
                if (p < 64) {
                    on = ((key.lo >> p) & 1ULL) != 0ULL;
                } else {
                    on = ((key.hi >> (p - 64)) & 1ULL) != 0ULL;
                }
                if (!on) {
                    continue;
                }

                const int i = pair_list[static_cast<std::size_t>(p)].first;
                const int j = pair_list[static_cast<std::size_t>(p)].second;
                adj[static_cast<std::size_t>(i)] |= static_cast<std::uint16_t>(1U << j);
                adj[static_cast<std::size_t>(j)] |= static_cast<std::uint16_t>(1U << i);
            }

            adj_masks.push_back(adj);
        }

        const int mask_count = 1 << n;
        std::vector<std::uint8_t> best(static_cast<std::size_t>(mask_count), 0);
        std::vector<std::uint8_t> score(static_cast<std::size_t>(mask_count), 0);

        for (const auto& adj : adj_masks) {
            score[0] = 0;
            for (int mask = 1; mask < mask_count; ++mask) {
                const int lb = mask & -mask;
                const int i = __builtin_ctz(static_cast<unsigned>(lb));
                const int prev = mask ^ lb;
                const std::uint16_t m = adj[static_cast<std::size_t>(i)];
                const int add = __builtin_popcount(static_cast<unsigned>(prev & m));
                score[static_cast<std::size_t>(mask)] =
                    static_cast<std::uint8_t>(score[static_cast<std::size_t>(prev)] + add);

                if (score[static_cast<std::size_t>(mask)] > best[static_cast<std::size_t>(mask)]) {
                    best[static_cast<std::size_t>(mask)] = score[static_cast<std::size_t>(mask)];
                }
            }
        }

        u64 numerator = 0;
        for (int mask = 0; mask < mask_count; ++mask) {
            const int consecutive_hh = __builtin_popcount(static_cast<unsigned>(mask & (mask >> 1)));
            numerator += static_cast<u64>(best[static_cast<std::size_t>(mask)]) +
                         static_cast<u64>(consecutive_hh);
        }
        return numerator;
    }
};

std::string solve_average_decimal(int n) {
    FoldingSolver solver(n);
    const u64 numerator = solver.solve_numerator();
    return FoldingSolver::exact_average_decimal(numerator, n);
}

bool run_checkpoints() {
    const std::string avg8 = solve_average_decimal(8);
    if (avg8 != "3.3203125") {
        std::cerr << "Checkpoint failed for n=8: got " << avg8 << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve_average_decimal(options.n) << '\n';
    return 0;
}

Python

def solve():
    n = 15

    grid_side = 2 * n + 1
    offset = n

    pair_list = []
    pair_index = {}
    idx = 0
    for i in range(n):
        for j in range(i + 2, n):
            pair_index[(i, j)] = idx
            pair_list.append((i, j))
            idx += 1

    num_pairs = len(pair_list)

    dx = [1, -1, 0, 0]
    dy = [0, 0, 1, -1]

    contact_maps = set()

    px = [0] * n
    py = [0] * n
    occupied = bytearray(grid_side * grid_side)

    def cell_id(x, y):
        return y * grid_side + x

    def add_contact_map():
        key = 0
        for p_idx in range(num_pairs):
            i, j = pair_list[p_idx]
            manhattan = abs(px[i] - px[j]) + abs(py[i] - py[j])
            if manhattan == 1:
                key |= (1 << p_idx)
        contact_maps.add(key)

    def dfs(idx_pos):
        if idx_pos == n - 1:
            add_contact_map()
            return
        cx = px[idx_pos]
        cy = py[idx_pos]
        for d in range(4):
            nx = cx + dx[d]
            ny = cy + dy[d]
            cid = cell_id(nx, ny)
            if occupied[cid]:
                continue
            occupied[cid] = 1
            px[idx_pos + 1] = nx
            py[idx_pos + 1] = ny
            dfs(idx_pos + 1)
            occupied[cid] = 0

    import sys
    sys.setrecursionlimit(500000)

    x0, y0 = offset, offset
    x1, y1 = offset + 1, offset
    occupied[cell_id(x0, y0)] = 1
    occupied[cell_id(x1, y1)] = 1
    px[0] = x0; py[0] = y0
    px[1] = x1; py[1] = y1
    dfs(1)

    maps_list = list(contact_maps)

    # Build adjacency masks
    adj_masks = []
    for key in maps_list:
        adj = [0] * n
        for p_idx in range(num_pairs):
            if (key >> p_idx) & 1:
                i, j = pair_list[p_idx]
                adj[i] |= (1 << j)
                adj[j] |= (1 << i)
        adj_masks.append(adj)

    mask_count = 1 << n
    best = bytearray(mask_count)

    popcount = [0] * mask_count
    for m in range(1, mask_count):
        popcount[m] = popcount[m >> 1] + (m & 1)

    for adj in adj_masks:
        score = bytearray(mask_count)
        for mask in range(1, mask_count):
            lb = mask & (-mask)
            i = lb.bit_length() - 1
            prev = mask ^ lb
            add = popcount[prev & adj[i]]
            score[mask] = score[prev] + add
            if score[mask] > best[mask]:
                best[mask] = score[mask]

    numerator = 0
    for mask in range(mask_count):
        consecutive_hh = popcount[mask & (mask >> 1)]
        numerator += best[mask] + consecutive_hh

    # exact_average_decimal
    denominator = 1 << n
    integer_part = numerator // denominator
    rem = numerator % denominator
    result = str(integer_part)
    if rem > 0:
        result += '.'
        while rem > 0:
            rem *= 10
            digit = rem // denominator
            result += str(digit)
            rem %= denominator
    return result

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler300 {
    static class ContactKey {
        long lo, hi;

        ContactKey(long lo, long hi) {
            this.lo = lo;
            this.hi = hi;
        }

        @Override
        public boolean equals(Object o) {
            ContactKey that = (ContactKey) o;
            return this.lo == that.lo && this.hi == that.hi;
        }

        @Override
        public int hashCode() {
            long x = lo * -7046045454385353131L; // 11400714819323198485ull
            long y = hi * -4417276706812531889L; // 14029467366897019727ull
            return (int) (x ^ (y >>> 1));
        }
    }

    static class FoldingSolver {
        int n;
        int gridSide;
        int offset;

        int[] px, py;
        byte[] occupied;

        int[] dx = { 1, -1, 0, 0 };
        int[] dy = { 0, 0, 1, -1 };

        List<int[]> pairList = new ArrayList<>();
        Set<ContactKey> contactMaps = new HashSet<>();

        FoldingSolver(int n) {
            this.n = n;
            gridSide = 2 * n + 1;
            offset = n;

            px = new int[n];
            py = new int[n];
            occupied = new byte[gridSide * gridSide];

            for (int i = 0; i < n; ++i) {
                for (int j = i + 2; j < n; ++j) {
                    pairList.add(new int[] { i, j });
                }
            }
        }

        int cellId(int x, int y) {
            return y * gridSide + x;
        }

        void addCurrentContactMap() {
            long lo = 0, hi = 0;
            for (int p = 0; p < pairList.size(); ++p) {
                int[] pair = pairList.get(p);
                int i = pair[0], j = pair[1];
                int manhattan = Math.abs(px[i] - px[j]) + Math.abs(py[i] - py[j]);
                if (manhattan == 1) {
                    if (p < 64)
                        lo |= (1L << p);
                    else
                        hi |= (1L << (p - 64));
                }
            }
            contactMaps.add(new ContactKey(lo, hi));
        }

        void dfs(int idx) {
            if (idx == n - 1) {
                addCurrentContactMap();
                return;
            }

            int cx = px[idx];
            int cy = py[idx];

            for (int dir = 0; dir < 4; ++dir) {
                int nx = cx + dx[dir];
                int ny = cy + dy[dir];
                int cid = cellId(nx, ny);
                if (occupied[cid] != 0)
                    continue;

                occupied[cid] = 1;
                px[idx + 1] = nx;
                py[idx + 1] = ny;
                dfs(idx + 1);
                occupied[cid] = 0;
            }
        }

        List<ContactKey> enumerateMaps() {
            int x0 = offset, y0 = offset;
            int x1 = offset + 1, y1 = offset;

            occupied[cellId(x0, y0)] = 1;
            occupied[cellId(x1, y1)] = 1;
            px[0] = x0;
            py[0] = y0;
            px[1] = x1;
            py[1] = y1;

            dfs(1);
            return new ArrayList<>(contactMaps);
        }

        long solveNumerator() {
            List<ContactKey> maps = enumerateMaps();
            List<int[]> adjMasks = new ArrayList<>(maps.size());

            for (ContactKey key : maps) {
                int[] adj = new int[n];
                for (int p = 0; p < pairList.size(); ++p) {
                    boolean on = false;
                    if (p < 64)
                        on = ((key.lo >>> p) & 1L) != 0;
                    else
                        on = ((key.hi >>> (p - 64)) & 1L) != 0;
                    if (on) {
                        int[] pair = pairList.get(p);
                        int i = pair[0], j = pair[1];
                        adj[i] |= (1 << j);
                        adj[j] |= (1 << i);
                    }
                }
                adjMasks.add(adj);
            }

            int maskCount = 1 << n;
            byte[] best = new byte[maskCount];
            byte[] score = new byte[maskCount];

            for (int[] adj : adjMasks) {
                score[0] = 0;
                for (int mask = 1; mask < maskCount; ++mask) {
                    int lb = mask & -mask;
                    int i = Integer.numberOfTrailingZeros(lb);
                    int prev = mask ^ lb;
                    int m = adj[i];
                    int add = Integer.bitCount(prev & m);
                    score[mask] = (byte) (score[prev] + add);
                    if (score[mask] > best[mask]) {
                        best[mask] = score[mask];
                    }
                }
            }

            long numerator = 0;
            for (int mask = 0; mask < maskCount; ++mask) {
                int consecutiveHh = Integer.bitCount(mask & (mask >>> 1));
                numerator += best[mask] + consecutiveHh;
            }
            return numerator;
        }
    }

    static String exactAverageDecimal(long numerator, int n) {
        long denominator = 1L << n;
        long integerPart = numerator / denominator;
        long rem = numerator % denominator;

        StringBuilder out = new StringBuilder();
        out.append(integerPart);
        if (rem == 0)
            return out.toString();

        out.append('.');
        while (rem != 0) {
            rem *= 10;
            long digit = rem / denominator;
            out.append(digit);
            rem %= denominator;
        }
        return out.toString();
    }

    public static String solve() {
        int n = 15;
        FoldingSolver solver = new FoldingSolver(n);
        long numerator = solver.solveNumerator();
        return exactAverageDecimal(numerator, n);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}