Problem 867: Tiling Dodecagon

View on Project Euler

Project Euler Problem 867 Solution

EulerSolve provides an optimized solution for Project Euler Problem 867, Tiling Dodecagon, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 867 asks for the number of admissible tilings of a dodecagon of order \(n\), and the implementations evaluate the case \(n=10\) modulo $$M=10^9+7.$$ The key point is that the geometry is not attacked tile by tile. Instead, the region is converted into a layered conflict graph, and a valid tiling becomes an independent set on that graph. Once that translation is made, the task becomes a transfer-matrix problem on rows whose lengths change by \(1\) at each step. Mathematical Approach For any finite sequence of row lengths \((\ell_1,\ell_2,\dots,\ell_t)\), let $$C(\ell_1,\ell_2,\dots,\ell_t)$$ denote the number of independent sets in the layered graph whose \(i\)-th row has \(\ell_i\) positions. The implementations solve the dodecagon problem by evaluating a small family of such row-profile counts and then combining them recursively. Step 1: Encode Each Row by an Independent Bitmask A row of length \(L\) is represented by a bitmask \(m\in\{0,1,\dots,2^L-1\}\). A bit \(1\) means that the corresponding position is selected in the independent set. Two adjacent positions in the same row cannot both be selected, so a legal mask must satisfy $$m \text{ is legal} \iff (m \mathbin{\&} (m \ll 1))=0.$$ Equivalently, a legal mask is an independent set of the path graph on \(L\) vertices....

Detailed mathematical approach

Problem Summary

Problem 867 asks for the number of admissible tilings of a dodecagon of order \(n\), and the implementations evaluate the case \(n=10\) modulo

$$M=10^9+7.$$

The key point is that the geometry is not attacked tile by tile. Instead, the region is converted into a layered conflict graph, and a valid tiling becomes an independent set on that graph. Once that translation is made, the task becomes a transfer-matrix problem on rows whose lengths change by \(1\) at each step.

Mathematical Approach

For any finite sequence of row lengths \((\ell_1,\ell_2,\dots,\ell_t)\), let

$$C(\ell_1,\ell_2,\dots,\ell_t)$$

denote the number of independent sets in the layered graph whose \(i\)-th row has \(\ell_i\) positions. The implementations solve the dodecagon problem by evaluating a small family of such row-profile counts and then combining them recursively.

Step 1: Encode Each Row by an Independent Bitmask

A row of length \(L\) is represented by a bitmask \(m\in\{0,1,\dots,2^L-1\}\). A bit \(1\) means that the corresponding position is selected in the independent set. Two adjacent positions in the same row cannot both be selected, so a legal mask must satisfy

$$m \text{ is legal} \iff (m \mathbin{\&} (m \ll 1))=0.$$

Equivalently, a legal mask is an independent set of the path graph on \(L\) vertices. The number of such masks is Fibonacci-like, which is why the method remains practical even when all masks up to width \(2n-1\) are precomputed.

Step 2: Describe Compatibility Between Consecutive Rows

Suppose the current row has length \(L_c\) and the next row has length \(L_n\), where the geometry guarantees \(L_n=L_c\pm1\). Fix a legal mask \(b\) on the next row. It forbids certain positions on the current row because any touching pair would violate independence.

If the next row is longer, \(L_n=L_c+1\), then a selected position in \(b\) blocks the position directly under it and the one immediately to its left. If the next row is shorter, \(L_n=L_c-1\), it blocks the position directly above it and the one immediately to its right. Writing \(\operatorname{Forb}(b)\) for that forbidden set, a previous-row mask \(a\) is compatible exactly when

$$a \subseteq \overline{\operatorname{Forb}(b)}.$$

Therefore, if \(\operatorname{dp}_{\mathrm{cur}}[a]\) counts all partial configurations ending with mask \(a\), then the next layer satisfies

$$\operatorname{dp}_{\mathrm{next}}[b]=\sum_{a\subseteq \overline{\operatorname{Forb}(b)}} \operatorname{dp}_{\mathrm{cur}}[a].$$

Step 3: Turn the Transition into a Subset-Sum Lookup

Computing the sum above separately for every \(b\) would repeat the same submask enumeration many times. The implementations avoid that by forming the subset zeta transform

$$Z(S)=\sum_{A\subseteq S}\operatorname{dp}_{\mathrm{cur}}[A].$$

Once \(Z\) has been computed for all subsets \(S\subseteq\{0,\dots,L_c-1\}\), every transition becomes a single table lookup:

$$\operatorname{dp}_{\mathrm{next}}[b]=Z\!\left(\overline{\operatorname{Forb}(b)}\right).$$

The standard bit-by-bit sweep computes \(Z\) in \(O(L_c2^{L_c})\) time, which is the dominant cost of each layer transition.

Step 4: Isolate the Two Canonical Profile Families

The full dodecagon is assembled from two recurring shapes, and both are counted by the same row-profile engine.

The first is a symmetric belt whose row lengths grow from \(n\) up to \(2n-1\) and then shrink back to \(n\):

$$B_n=C(n,n+1,\dots,2n-1,2n-2,\dots,n).$$

The second is a corner staircase of height \(h-1\):

$$K_{n,h}=C(n-2,n-3,\dots,n-h).$$

So the geometric complexity of the dodecagon is reduced to repeatedly evaluating these belt and corner profile counts.

Step 5: Reassemble the Dodecagon Recursively

Let \(Q(u,v)\) be the memoized count for the reduced central configuration that remains after stripping equal corner layers. The implementations use the recurrence

$$Q(u,0)=B_u,$$

$$Q(u,v)=\mathbf{1}_{(u,v)=(1,1)}+\sum_{w=0}^{u-1} Q(v,w)\,K_{u,u-w}^{\,6} \qquad (v>0).$$

The exponent \(6\) comes from the six congruent corner sectors of the dodecagon. After the interface parameter \(w\) is fixed, those six sectors contribute independently, so their counts multiply. The target quantity is then recovered by

$$A_n=2Q(n,n)-\mathbf{1}_{n=1}.$$

This formula is exactly the one implemented in all three languages.

Worked Example: The Cases \(n=1\) and \(n=2\)

For \(n=1\), the symmetric belt consists of a single row of length \(1\), so \(B_1=C(1)=2\): the mask may be empty or occupied. The corner profile is empty, hence \(K_{1,1}=1\). Therefore

$$Q(1,1)=1+Q(1,0)\,K_{1,1}^{\,6}=1+2=3,$$

and the final count is

$$A_1=2\cdot 3-1=5.$$

For \(n=2\), the belt profile is \((2,3,2)\), and the implementation evaluates

$$B_2=C(2,3,2)=19.$$

Both relevant corner staircases contribute \(1\), so

$$Q(2,1)=Q(1,0)+Q(1,1)=2+3=5,$$

$$Q(2,2)=Q(2,0)+Q(2,1)=19+5=24,$$

and finally

$$A_2=2\cdot 24=48.$$

These two values are the small checkpoints used by the implementations, so they provide a good sanity check for the whole derivation.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they precompute every legal row mask for each width from \(0\) to \(2n-1\). Next they run the row-profile dynamic program described above, caching results by the row-length sequence so that repeated belt or corner profiles are evaluated only once. On each transition they build the subset-sum table for the current row, use it to fill the next row in one pass over the legal masks, and then sum the counts on the last row.

Above that transfer-matrix layer, the implementation memoizes the one-parameter belt counts, the two-parameter corner counts, and the two-parameter recursive assembly values. Modular exponentiation is used only when the six identical corner contributions must be raised to the sixth power. Because the same subproblems occur many times, memoization is essential to keeping the recursive stage small.

Complexity Analysis

Let \(W=2n-1\) be the maximum row width. Precomputing all legal masks costs \(O(2^W)\) time and storage. For a single row-profile sequence \((\ell_1,\dots,\ell_t)\), the transfer stage costs

$$O\!\left(\sum_{i=1}^{t-1} \ell_i 2^{\ell_i}\right),$$

because each layer transition is dominated by one subset-zeta transform on a width-\(\ell_i\) row. The recursive assembly has \(O(n^2)\) memo states, and each state sums over at most \(u\) interface values, so the extra gluing work is \(O(n^3)\) after cached profile counts are available. In practice the profile DP dominates, while the memory footprint is \(O(2^W)\) for the active DP arrays plus the memo tables.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=867
  2. Independent set: Wikipedia - Independent set (graph theory)
  3. Transfer-matrix method: Wikipedia - Transfer-matrix method
  4. Dynamic programming: Wikipedia - Dynamic programming
  5. Subset enumeration and subset DP: cp-algorithms - Enumerating submasks of a bitmask

Problem 867 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <functional>
#include <iostream>
#include <map>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

using i64 = std::int64_t;
using u64 = std::uint64_t;

static constexpr i64 MOD = 1'000'000'007LL;

static i64 mod_pow(i64 a, i64 e) {
    i64 r = 1;
    while (e > 0) {
        if (e & 1LL) r = static_cast<i64>((static_cast<__int128>(r) * a) % MOD);
        a = static_cast<i64>((static_cast<__int128>(a) * a) % MOD);
        e >>= 1LL;
    }
    return r;
}

class Solver867 {
public:
    explicit Solver867(int n) : n_(n), max_len_(2 * n - 1), indep_(max_len_ + 1) {
        build_independent_masks();
    }

    i64 T(int n) {
        i64 ans = (2LL * R(n, n)) % MOD;
        if (n == 1) ans = (ans + MOD - 1) % MOD;
        return ans;
    }

private:
    int n_;
    int max_len_;
    std::vector<std::vector<int>> indep_;
    std::map<std::string, i64> count_cache_;
    std::unordered_map<int, i64> h_cache_;
    std::unordered_map<u64, i64> f_cache_;
    std::unordered_map<u64, i64> r_cache_;

    static u64 pair_key(int a, int b) {
        return (static_cast<u64>(static_cast<std::uint32_t>(a)) << 32) |
               static_cast<u64>(static_cast<std::uint32_t>(b));
    }

    void build_independent_masks() {
        for (int L = 0; L <= max_len_; ++L) {
            int lim = 1 << L;
            indep_[L].reserve(lim);
            for (int m = 0; m < lim; ++m) {
                if ((m & (m << 1)) == 0) indep_[L].push_back(m);
            }
        }
    }

    std::vector<i64> zeta_subset_sum(const std::vector<i64>& v, int L) const {
        std::vector<i64> out = v;
        int size = 1 << L;
        for (int bit = 0; bit < L; ++bit) {
            int step = 1 << bit;
            int block = step << 1;
            for (int start = 0; start < size; start += block) {
                int mid = start + step;
                int end = start + block;
                for (int m = mid; m < end; ++m) {
                    i64 x = out[m] + out[m - step];
                    out[m] = (x >= MOD ? x - MOD : x);
                }
            }
        }
        return out;
    }

    static std::string encode_lengths(const std::vector<int>& lens) {
        std::string key;
        key.reserve(lens.size() * 2 + 1);
        for (int x : lens) {
            key.push_back(static_cast<char>(x + 1));
            key.push_back('|');
        }
        return key;
    }

    i64 count_independent_sets(const std::vector<int>& row_lengths) {
        std::string key = encode_lengths(row_lengths);
        auto it = count_cache_.find(key);
        if (it != count_cache_.end()) return it->second;

        if (row_lengths.empty()) {
            count_cache_[key] = 1;
            return 1;
        }

        int L0 = row_lengths[0];
        std::vector<i64> dp(1 << L0, 0);
        for (int m : indep_[L0]) dp[m] = 1;

        for (std::size_t i = 1; i < row_lengths.size(); ++i) {
            int Lc = row_lengths[i - 1];
            int Ln = row_lengths[i];

            std::vector<i64> subs = zeta_subset_sum(dp, Lc);
            std::vector<i64> ndp(1 << Ln, 0);
            int full = (1 << Lc) - 1;

            if (Ln == Lc + 1) {
                for (int b : indep_[Ln]) {
                    int forb = (b | (b >> 1)) & full;
                    int allowed = full ^ forb;
                    ndp[b] = subs[allowed];
                }
            } else {
                for (int b : indep_[Ln]) {
                    int forb = (b | (b << 1)) & full;
                    int allowed = full ^ forb;
                    ndp[b] = subs[allowed];
                }
            }

            dp.swap(ndp);
        }

        int lastL = row_lengths.back();
        i64 total = 0;
        for (int m : indep_[lastL]) {
            total += dp[m];
            if (total >= MOD) total -= MOD;
        }

        count_cache_[key] = total;
        return total;
    }

    i64 H(int n) {
        auto it = h_cache_.find(n);
        if (it != h_cache_.end()) return it->second;

        std::vector<int> lens;
        lens.reserve(2 * n - 1);
        for (int x = n; x <= 2 * n - 1; ++x) lens.push_back(x);
        for (int x = 2 * n - 2; x >= n; --x) lens.push_back(x);

        i64 val = count_independent_sets(lens);
        h_cache_[n] = val;
        return val;
    }

    i64 F(int n, int h) {
        u64 key = pair_key(n, h);
        auto it = f_cache_.find(key);
        if (it != f_cache_.end()) return it->second;

        int rows = h - 1;
        std::vector<int> lens;
        lens.reserve(std::max(0, rows));
        for (int i = 0; i < rows; ++i) {
            lens.push_back(std::max(0, (n - 2) - i));
        }

        i64 val = count_independent_sets(lens);
        f_cache_[key] = val;
        return val;
    }

    i64 R(int u, int v) {
        u64 key = pair_key(u, v);
        auto it = r_cache_.find(key);
        if (it != r_cache_.end()) return it->second;

        i64 res;
        if (v == 0) {
            res = H(u);
        } else {
            res = (u == 1 && v == 1) ? 1 : 0;
            for (int w = 0; w < u; ++w) {
                i64 corner = F(u, u - w);
                i64 add = static_cast<i64>((static_cast<__int128>(R(v, w)) * mod_pow(corner, 6)) % MOD);
                res += add;
                if (res >= MOD) res -= MOD;
            }
        }

        r_cache_[key] = res;
        return res;
    }
};

int main() {
    Solver867 solver(10);

    assert(solver.T(1) == 5);
    assert(solver.T(2) == 48);

    std::cout << solver.T(10) << '\n';
    return 0;
}

Python

MOD = 1000000007

def mod_pow(a, e):
    return pow(a, e, MOD)

class Solver867:
    def __init__(self, n):
        self.n = n
        self.max_len = 2 * n - 1
        self.indep = [[] for _ in range(self.max_len + 1)]
        self.build_independent_masks()
        self.count_cache = {}
        self.h_cache = {}
        self.f_cache = {}
        self.r_cache = {}

    def pair_key(self, a, b):
        return (a << 32) | b

    def build_independent_masks(self):
        for L in range(self.max_len + 1):
            lim = 1 << L
            for m in range(lim):
                if (m & (m << 1)) == 0:
                    self.indep[L].append(m)

    def zeta_subset_sum(self, v, L):
        out = list(v)
        size = 1 << L
        for bit in range(L):
            step = 1 << bit
            block = step << 1
            for start in range(0, size, block):
                mid = start + step
                end = start + block
                for m in range(mid, end):
                    x = out[m] + out[m - step]
                    if x >= MOD:
                        x -= MOD
                    out[m] = x
        return out

    def encode_lengths(self, lens):
        return "|".join(map(str, lens))

    def count_independent_sets(self, row_lengths):
        key = self.encode_lengths(row_lengths)
        if key in self.count_cache:
            return self.count_cache[key]

        if not row_lengths:
            self.count_cache[key] = 1
            return 1

        L0 = row_lengths[0]
        dp = [0] * (1 << L0)
        for m in self.indep[L0]:
            dp[m] = 1

        for i in range(1, len(row_lengths)):
            Lc = row_lengths[i - 1]
            Ln = row_lengths[i]

            subs = self.zeta_subset_sum(dp, Lc)
            ndp = [0] * (1 << Ln)
            full = (1 << Lc) - 1

            if Ln == Lc + 1:
                for b in self.indep[Ln]:
                    forb = (b | (b >> 1)) & full
                    allowed = full ^ forb
                    ndp[b] = subs[allowed]
            else:
                for b in self.indep[Ln]:
                    forb = (b | (b << 1)) & full
                    allowed = full ^ forb
                    ndp[b] = subs[allowed]

            dp = ndp

        lastL = row_lengths[-1]
        total = 0
        for m in self.indep[lastL]:
            total += dp[m]
            if total >= MOD:
                total -= MOD

        self.count_cache[key] = total
        return total

    def H(self, n):
        if n in self.h_cache:
            return self.h_cache[n]

        lens = []
        for x in range(n, 2 * n):
            lens.append(x)
        for x in range(2 * n - 2, n - 1, -1):
            lens.append(x)

        val = self.count_independent_sets(lens)
        self.h_cache[n] = val
        return val

    def F(self, n, h):
        key = self.pair_key(n, h)
        if key in self.f_cache:
            return self.f_cache[key]

        rows = h - 1
        lens = []
        for i in range(max(0, rows)):
            lens.append(max(0, (n - 2) - i))

        val = self.count_independent_sets(lens)
        self.f_cache[key] = val
        return val

    def R(self, u, v):
        key = self.pair_key(u, v)
        if key in self.r_cache:
            return self.r_cache[key]

        if v == 0:
            res = self.H(u)
        else:
            res = 1 if (u == 1 and v == 1) else 0
            for w in range(u):
                corner = self.F(u, u - w)
                add = (self.R(v, w) * mod_pow(corner, 6)) % MOD
                res += add
                if res >= MOD:
                    res -= MOD

        self.r_cache[key] = res
        return res

    def T(self, n):
        ans = (2 * self.R(n, n)) % MOD
        if n == 1:
            ans = (ans + MOD - 1) % MOD
        return ans

def solve():
    solver = Solver867(10)
    return str(solver.T(10))

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

Java

import java.util.ArrayList;
import java.util.HashMap;

public class Euler867 {
    static final long MOD = 1000000007L;

    static long modPow(long a, long e) {
        long r = 1;
        while (e > 0) {
            if ((e & 1) == 1)
                r = (r * a) % MOD;
            a = (a * a) % MOD;
            e >>= 1;
        }
        return r;
    }

    static class Solver867 {
        int n;
        int maxLen;
        ArrayList<ArrayList<Integer>> indep;
        HashMap<String, Long> countCache;
        HashMap<Integer, Long> hCache;
        HashMap<Long, Long> fCache;
        HashMap<Long, Long> rCache;

        Solver867(int n) {
            this.n = n;
            this.maxLen = 2 * n - 1;
            this.indep = new ArrayList<>();
            for (int i = 0; i <= maxLen; ++i) {
                indep.add(new ArrayList<>());
            }
            this.countCache = new HashMap<>();
            this.hCache = new HashMap<>();
            this.fCache = new HashMap<>();
            this.rCache = new HashMap<>();
            buildIndependentMasks();
        }

        long pairKey(int a, int b) {
            return (((long) a) << 32) | ((long) b & 0xFFFFFFFFL);
        }

        void buildIndependentMasks() {
            for (int L = 0; L <= maxLen; ++L) {
                int lim = 1 << L;
                for (int m = 0; m < lim; ++m) {
                    if ((m & (m << 1)) == 0)
                        indep.get(L).add(m);
                }
            }
        }

        long[] zetaSubsetSum(long[] v, int L) {
            long[] out = v.clone();
            int size = 1 << L;
            for (int bit = 0; bit < L; ++bit) {
                int step = 1 << bit;
                int block = step << 1;
                for (int start = 0; start < size; start += block) {
                    int mid = start + step;
                    int end = start + block;
                    for (int m = mid; m < end; ++m) {
                        long x = out[m] + out[m - step];
                        out[m] = (x >= MOD ? x - MOD : x);
                    }
                }
            }
            return out;
        }

        String encodeLengths(ArrayList<Integer> lens) {
            StringBuilder sb = new StringBuilder();
            for (int x : lens) {
                sb.append(x).append("|");
            }
            return sb.toString();
        }

        long countIndependentSets(ArrayList<Integer> rowLengths) {
            String key = encodeLengths(rowLengths);
            if (countCache.containsKey(key))
                return countCache.get(key);

            if (rowLengths.isEmpty()) {
                countCache.put(key, 1L);
                return 1L;
            }

            int L0 = rowLengths.get(0);
            long[] dp = new long[1 << L0];
            for (int m : indep.get(L0))
                dp[m] = 1;

            for (int i = 1; i < rowLengths.size(); ++i) {
                int Lc = rowLengths.get(i - 1);
                int Ln = rowLengths.get(i);

                long[] subs = zetaSubsetSum(dp, Lc);
                long[] ndp = new long[1 << Ln];
                int full = (1 << Lc) - 1;

                if (Ln == Lc + 1) {
                    for (int b : indep.get(Ln)) {
                        int forb = (b | (b >> 1)) & full;
                        int allowed = full ^ forb;
                        ndp[b] = subs[allowed];
                    }
                } else {
                    for (int b : indep.get(Ln)) {
                        int forb = (b | (b << 1)) & full;
                        int allowed = full ^ forb;
                        ndp[b] = subs[allowed];
                    }
                }

                dp = ndp;
            }

            int lastL = rowLengths.get(rowLengths.size() - 1);
            long total = 0;
            for (int m : indep.get(lastL)) {
                total += dp[m];
                if (total >= MOD)
                    total -= MOD;
            }

            countCache.put(key, total);
            return total;
        }

        long H(int n) {
            if (hCache.containsKey(n))
                return hCache.get(n);

            ArrayList<Integer> lens = new ArrayList<>();
            for (int x = n; x <= 2 * n - 1; ++x)
                lens.add(x);
            for (int x = 2 * n - 2; x >= n; --x)
                lens.add(x);

            long val = countIndependentSets(lens);
            hCache.put(n, val);
            return val;
        }

        long F(int n, int h) {
            long key = pairKey(n, h);
            if (fCache.containsKey(key))
                return fCache.get(key);

            int rows = h - 1;
            ArrayList<Integer> lens = new ArrayList<>();
            for (int i = 0; i < Math.max(0, rows); ++i) {
                lens.add(Math.max(0, (n - 2) - i));
            }

            long val = countIndependentSets(lens);
            fCache.put(key, val);
            return val;
        }

        long R(int u, int v) {
            long key = pairKey(u, v);
            if (rCache.containsKey(key))
                return rCache.get(key);

            long res;
            if (v == 0) {
                res = H(u);
            } else {
                res = (u == 1 && v == 1) ? 1 : 0;
                for (int w = 0; w < u; ++w) {
                    long corner = F(u, u - w);
                    long add = (R(v, w) * modPow(corner, 6)) % MOD;
                    res += add;
                    if (res >= MOD)
                        res -= MOD;
                }
            }

            rCache.put(key, res);
            return res;
        }

        long T(int n) {
            long ans = (2L * R(n, n)) % MOD;
            if (n == 1)
                ans = (ans + MOD - 1) % MOD;
            return ans;
        }
    }

    public static String solve() {
        Solver867 solver = new Solver867(10);
        return Long.toString(solver.T(10));
    }

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