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

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.

For example, the solver finds only

$$41$$

distinct contact maps for $n=8$, and

$$12495$$

distinct contact maps for $n=15$.

3. Why consecutive H-H pairs are added separately

If positions $i$ and $i+1$ are both H, then they are always adjacent in the chain, no matter how the fold bends. Therefore that contribution is independent of the contact map. For an HP mask $h\in\{0,1\}^n$, the number of consecutive H-H bonds is

$$\operatorname{popcount}(h\ \&\ (h\gg 1)).$$

This explains the extra term in the code:

$$\text{total}(h,C)=\text{nonconsecutive\_contacts}(h,C)+\operatorname{popcount}(h\ \&\ (h\gg 1)).$$

So the optimization over folds only needs to maximize the non-consecutive part.

4. Scoring one contact map by subset DP

Fix one contact map $C$. The code converts it to adjacency bitmasks $N_C(i)$, where bit $j$ is set if $(i,j)$ is a non-consecutive contact in that fold. For an H-mask $m$, define $\text{score}_C(m)$ as the number of contact edges whose two endpoints are both H.

Instead of checking every contact edge from scratch for every mask, the implementation uses the recurrence

$$\text{score}_C(m)=\text{score}_C(m\setminus\{i\})+\bigl|N_C(i)\cap(m\setminus\{i\})\bigr|,$$

where $i$ is the least significant set bit of $m$. This works because when $i$ is inserted last, exactly the contacts from $i$ to already-present H positions are new. Every contact is counted once and only once.

5. Best fold for each HP sequence

Let $B(m)$ be the best non-consecutive contact count over all unique maps:

$$B(m)=\max_C \text{score}_C(m).$$

The code maintains an array best[mask] and updates it while scanning all deduplicated contact maps. After that, the full optimal score for sequence mask $m$ is

$$B(m)+\operatorname{popcount}(m\ \&\ (m\gg 1)).$$

6. Exact averaging over all $2^n$ sequences

The final numerator is

$$\sum_{m=0}^{2^n-1}\left(B(m)+\operatorname{popcount}(m\ \&\ (m\gg 1))\right).$$

Since every HP sequence is equally likely, the desired expectation is

$$\frac{1}{2^n}\sum_{m=0}^{2^n-1}\left(B(m)+\operatorname{popcount}(m\ \&\ (m\gg 1))\right).$$

The denominator is a power of two, so the result is converted to an exact terminating binary-rational decimal by long division.

7. Checkpoints

The published checkpoint is

$$n=8 \quad \Longrightarrow \quad \frac{850}{256}=3.3203125.$$

The full solver for $n=15$ returns

$$8.0540771484375.$$

How the Code Works

dfs(...) enumerates all self-avoiding walks after fixing the first edge. add_current_contact_map() converts one finished walk to a packed bit key and inserts it into a set. Then each unique map is turned into adjacency masks, and a subset DP computes $\text{score}_C(m)$ for all $m$. The array best stores the maximum over all maps, after which the sequence-independent consecutive-HH term is added and averaged.

Complexity Analysis

If $M$ is the number of unique contact maps, the scoring phase costs

$$O(M\cdot 2^n)$$

with very small constants because every update is just a few bit operations. Memory is

$$O(M+2^n).$$

The expensive combinatorial step is the self-avoiding walk enumeration, but for $n=15$ it is still manageable once symmetry is fixed and maps are deduplicated.

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());
    }
}