Problem 364: Comfortable Distance

Project Euler Problem 364 Solution

EulerSolve provides an optimized solution for Project Euler Problem 364, Comfortable Distance, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary There are $n$ seats in one row, and people sit down one after another. At every step, the allowed choices are exactly the empty seats with the smallest possible number of occupied neighbours. If $T(n)$ denotes the number of complete seating orders, the program must compute $T(10^6) \bmod 100000007$. A direct state search is exponential, so the solution counts orders through the gap pattern that appears after the last seat with zero occupied neighbours has been taken. Mathematical Approach Step 1: The greedy rule forces three phases The neighbour count of an empty seat can only be $0$, $1$, or $2$. Because the rule always chooses the minimum, every valid seating order is split into three consecutive phases: first all choices with $0$ occupied neighbours, then all choices with $1$ occupied neighbour, and finally all choices with $2$ occupied neighbours. No move from a later phase can occur while an earlier phase is still available. Step 2: The zero-neighbour skeleton Stop the process exactly when no seat with $0$ occupied neighbours remains. The occupied seats are then isolated singletons; let their number be $m$. Between consecutive singletons there are $m-1$ internal gaps. Each such gap must have length $1$ or $2$, because a longer internal gap would still contain a seat with zero occupied neighbours....

Detailed mathematical approach

Problem Summary

A row has $N$ seats, numbered $1$ to $N$, and $N$ people enter one after another. Each newcomer applies the first applicable rule: take a seat whose adjacent seat or seats are all unoccupied; otherwise take a seat for which exactly one adjacent seat is occupied; otherwise take any remaining seat. $T(N)$ is the number of possibilities, that is of sequences of choices that fill all $N$ seats. The statement's figure shows the $T(4)=8$ possibilities, and its checkpoints are $T(10)=61632$ together with

$$T(1000)\equiv47255094\pmod{100000007}.$$

We must find $T(10^6)\bmod 100000007$. ^[1]

The whole solution in one sentence: every history breaks into three phases — first the seats that have no occupied neighbour, then those with one, then those with two — so a history is described by the set of seats filled in the first phase plus an order inside each phase, and counting those sets by their pattern of gaps turns $T(N)$ into a sum of about $N/2$ products of factorials.

Mathematical Approach

1. The rules say "fewest occupied neighbours"

For an empty seat $i$ let $d(i)$ be the number of its neighbours that are occupied. For $N\ge2$ the two end seats have one neighbour and every other seat has two, so $d(i)\in\{0,1,2\}$. Rule 1 applies when some empty seat has $d=0$, rule 2 when the smallest value of $d$ is $1$, and rule 3 otherwise. The three rules together therefore say one thing: the newcomer takes any empty seat at which $d$ is smallest.

Three seats already show that this is a real restriction. At the start every seat is free of occupied neighbours, so the first person has three choices. After seat $2$, seats $1$ and $3$ both have $d=1$, the second person may take either and the third is then forced: $2$ possibilities. After seat $1$, seat $3$ still has $d=0$ while seat $2$ has $d=1$, so rule 1 forces seat $3$ next and seat $2$ last: $1$ possibility, and likewise after seat $3$. Hence $T(3)=4$, and the natural-looking order $1,2,3$ is not allowed.

2. The move types come in a fixed order

Call a move a $0$-move, a $1$-move or a $2$-move according to the value of $d$ at the seat that is taken; by Section 1 that value is the minimum of $d$ over the empty seats just before the move.

Claim. Consecutive move types never decrease.

Proof. Let $E$ and $E'$ be the empty seats before and after a move, so $E'\subseteq E$, and let $d'$ be the new counts; filling a seat cannot lower $d$ at a seat that stays empty, so $d'\ge d$. While $E'$ is not empty, its minimum of $d'$ is at least its minimum of $d$, which is at least the minimum of $d$ over $E$; after the last move no next type exists. $\square$

So the $N$ moves split into phase 1 (the $0$-moves), phase 2 (the $1$-moves) and phase 3 (the $2$-moves), in that order. Phase 1 is never empty, because the first person faces an empty row, but the other two phases may be: for $N=3$ the history $1,3,2$ has types $0,0,2$ and no phase 2, while $2,1,3$ has types $0,1,1$ and no phase 3. An empty phase admits exactly one order, which the factorials below count as $0!=1$.

3. Phase 1 ends at a maximal independent set

Let $S$ be the set of seats taken during phase 1 and $m=|S|$.

No two seats of $S$ are adjacent: of two adjacent seats, the one taken later would have had $d\ge1$, so its move would not be a $0$-move.
Every seat outside $S$ has a neighbour in $S$: phase 1 ends only when no empty seat has $d=0$, and at that moment the occupied seats are exactly $S$.

In graph language the seats and the neighbouring pairs form a path $P_N$, and $S$ is an independent set that is also dominating, that is a maximal independent set: no seat can be added without creating an adjacent pair. ^[2] We call it the skeleton of the history, and since every history has exactly one skeleton, we may count histories skeleton by skeleton.

4. The arithmetic of the gaps

Write the skeleton as $p_1\lt p_2\lt\cdots\lt p_m$ and look at the blocks of empty seats it leaves: the left gap $\ell=p_1-1$, the right gap $r=N-p_m$, and the $m-1$ internal gaps $c_i=p_{i+1}-p_i-1$. Independence gives $c_i\ge1$; domination bounds them above, since an internal gap of three or more empty seats would contain a seat with two empty neighbours, and a boundary gap of two or more would leave the end seat with its only neighbour empty. Hence

$$c_i\in\{1,2\},\qquad \ell,r\in\{0,1\}.$$

Let $g_2$ be the number of internal gaps of length $2$ and put $s=\ell+r\in\{0,1,2\}$. Counting all seats once gives

$$N=m+\bigl((m-1)+g_2\bigr)+s=2m-1+s+g_2,$$

so $g_2$ is not free once $m$ and $s$ are chosen:

$$g_2=N-(2m-1+s),\qquad 0\le g_2\le m-1.$$

5. How many skeletons a pair $(m,s)$ carries

For $s=0$ the two boundary gaps are forced to be empty, for $s=2$ both have length $1$, and for $s=1$ there are two patterns, $(\ell,r)=(1,0)$ and $(0,1)$; they give different skeletons because $p_1$ differs. The multiplicities are therefore

$$w_0=1,\qquad w_1=2,\qquad w_2=1.$$

Inside, the only remaining freedom is which of the $m-1$ internal gaps are the long ones, so the number of skeletons with parameters $(m,s)$ is $w_s\binom{m-1}{g_2}$, a binomial coefficient: the number of ways to choose $g_2$ items out of $m-1$.

Summing this over $m$ and $s$ must give the total number of maximal independent sets of the path, and it does: for $N=1,\dots,10$ the totals are $1, 2, 2, 3, 4, 5, 7, 9, 12, 16$, the Padovan numbers, which satisfy $a(N)=a(N-2)+a(N-3)$. ^[2] ^[3] They grow like $\rho^N$ with the plastic number $\rho=1.3247179572\ldots$, the real root of $\rho^3=\rho+1$, ^[3] so at $N=10^6$ the skeletons cannot be listed one by one, and the pair $(m,s)$ is what lets us count them in bulk.

6. The orders inside phase 1

Fix a skeleton $S$. The histories whose phase-1 set is $S$ are exactly those whose first $m$ moves take the seats of $S$ in some order, and all $m!$ orders occur: no seat of $S$ is adjacent to another, so when one is taken its neighbours are still empty and the move is a legal $0$-move, and while some seat of $S$ is still empty it has $d=0$, so the minimum is $0$ and no later phase can begin.

7. The orders inside phase 2

At the end of phase 1 the empty seats are precisely the gaps, and their values of $d$ are easy to read off:

an internal gap of length $1$: its single seat has two occupied neighbours, $d=2$;
an internal gap of length $2$: each of its two seats has one occupied and one empty neighbour, $d=1$;
a boundary gap of length $1$: its end seat has its only neighbour occupied, $d=1$.

So phase 2 starts with $2g_2+s$ available seats. Taking one seat of a long internal gap leaves its partner with two occupied neighbours, hence with $d=2$, so no move ever creates a new seat with $d=1$; and different gaps never influence each other, because occupied seats separate them. Each long internal gap therefore contributes exactly one phase-2 move, with two choices, and each boundary gap of length $1$ one forced phase-2 move. The $s+g_2$ moves may be interleaved in any order, which gives

$$(s+g_2)!\,2^{g_2}$$

possibilities for the whole of phase 2.

8. The orders inside phase 3

After phase 2 every internal gap holds exactly one empty seat, with both neighbours occupied, and the boundary gaps are gone. These $m-1$ seats are pairwise non-adjacent, so filling one does not change the others, and every one of the $(m-1)!$ orders is a valid phase 3.

9. The formula

Multiplying the number of skeletons by the orders of the three phases gives the contribution of one pair $(m,s)$,

$$C(m,s)=w_s\binom{m-1}{g_2}m!\,(m-1)!\,(s+g_2)!\,2^{g_2},$$

with $g_2=N-(2m-1+s)$, and

$$T(N)=\sum_{m\ge1}\ \sum_{s=0}^{2}C(m,s),$$

where only the pairs with $0\le g_2\le m-1$ are kept. Every history has exactly one skeleton by Section 3, every skeleton really occurs by Section 6, and the three phase counts depend only on $m$, $s$ and $g_2$, so these contributions partition the histories and nothing is counted twice. For each $s$ the surviving values of $m$ form an interval:

$$\begin{aligned}m&\ge\lceil (N+2-s)/3\rceil,\\ m&\le\lfloor (N+1-s)/2\rfloor.\end{aligned}$$

At $N=10^6$ exactly $500\,000$ of the $3\,000\,000$ candidate pairs survive, with $m$ between $333\,334$ and $500\,000$.

The smallest cases come out right. For $N=3$ the pairs $(1,2)$ and $(2,0)$ contribute $2$ each, the $T(3)=4$ of Section 1. For $N=4$ the pair $(2,0)$ has $g_2=1$ and contributes $1\cdot1\cdot2!\cdot1!\cdot1!\cdot2=4$, the pair $(2,1)$ has $g_2=0$ and contributes $2\cdot1\cdot2!\cdot1!\cdot1!\cdot1=4$, and $T(4)=8$, the value of the figure.

10. A worked example: $T(10)=61632$

For $N=10$ five pairs survive. Each line lists $(m,s)$, then $g_2$, the number of skeletons $w_s\binom{m-1}{g_2}$, the orders $m!\,(s+g_2)!\,2^{g_2}(m-1)!$ of one skeleton, and the product:

$(4,0)$: $g_2=3$, $1$ skeleton, $24\cdot6\cdot8\cdot6=6912$ orders, $6912$;
$(4,1)$: $g_2=2$, $6$ skeletons, $24\cdot6\cdot4\cdot6=3456$ orders, $20736$;
$(4,2)$: $g_2=1$, $3$ skeletons, $24\cdot6\cdot2\cdot6=1728$ orders, $5184$;
$(5,0)$: $g_2=1$, $4$ skeletons, $120\cdot1\cdot2\cdot24=5760$ orders, $23040$;
$(5,1)$: $g_2=0$, $2$ skeletons, $120\cdot1\cdot1\cdot24=2880$ orders, $5760$.

The skeletons number $1+6+3+4+2=16$, the Padovan value of Section 5, and the products sum to the checkpoint of the statement:

$$\begin{aligned}&6912+20736+5184\\ &{}+23040+5760=61632.\end{aligned}$$

One of the six skeletons of the line $(4,1)$ is $\{2,5,8,10\}$: here $\ell=1$, $r=0$, and the internal gaps have lengths $2,2,1$. Phase 1 fills seats $2,5,8,10$ in any of $4!=24$ orders. Phase 2 has three moves: seat $1$, one of $\{3,4\}$, one of $\{6,7\}$; seat $9$ is excluded because it already has two occupied neighbours. That is $3!\cdot2^2=24$ possibilities. Phase 3 fills the three seats that are left, one per internal gap, in $3!=6$ orders. The skeleton carries $24\cdot24\cdot6=3456$ histories, and the six skeletons of that line give $20736$, as listed above.

11. Working modulo $100\,000\,007$

The modulus $P=100\,000\,007=10^8+7$ is prime, as trial division by every integer up to $10^4$ confirms. Since $P$ is far larger than $N+1$, no factorial of an index at most $N+1$ is divisible by $P$, so every factorial we need is invertible modulo $P$: an inverse of $a$ is a residue $b$ with $ab\equiv1\pmod P$, and multiplying by $b$ is how one divides by $a$. Fermat's little theorem states that $a^{P-1}\equiv1\pmod P$ whenever $a$ is not divisible by $P$, ^[6] so $a^{P-2}$ is the inverse of $a$. One exponentiation by squaring, which uses $O(\log P)$ modular multiplications, ^[7] gives the inverse of $(N+1)!$, and then the identity

$$(k!)^{-1}\equiv(k+1)\cdot\bigl((k+1)!\bigr)^{-1}$$

walks down through all inverse factorials in one pass. Since $\binom{m-1}{g_2}$ is $(m-1)!$ divided by $g_2!$ and by $(m-1-g_2)!$, its residue is one factorial times two inverse factorials, so a binomial coefficient costs two modular multiplications once the tables are ready. ^[8] A pleasant side effect of this unusually small modulus is that $P\lt2^{27}$, so a product of two residues stays below $2^{54}$ and fits into signed 64-bit arithmetic.

How the Code Works

mod_pow is exponentiation by squaring modulo $P$: it squares the base and multiplies it into the result whenever the lowest bit of the exponent is set. ^[7] C++ casts every product to __int128, Java keeps long arithmetic, safe by Section 11, and the Python wrapper calls the built-in three-argument pow.

comfortable_arrangements_mod(n), which Java spells comfortableArrangementsMod like all its identifiers, does the work. It fills fact with the factorials up to $n+1$, seeds inv_fact at index $n+1$ with mod_pow(fact[n+1], kMod - 2), runs the downward recurrence of Section 11, and tabulates the powers of two in pow2. The arrays boundary_sum and boundary_mult hold $(0,1,2)$ and $(1,2,1)$: the $s$ and the $w_s$ of Section 5.

The main loop runs over $m=1,\dots,n$ and, inside, the three boundary cases. For each pair it computes g2 as $n-(2m-1+s)$ and skips the pair unless $g_2$ lies in the interval of Section 4; at $n=10^6$ it visits all $3\,000\,000$ candidate pairs and evaluates the $500\,000$ surviving terms of Section 9. The coefficient $\binom{m-1}{g_2}$ is formed as fact[m-1] times inv_fact[g2] times inv_fact[m-1-g2], and the multiplicity, fact[m], fact[m-1], fact[s+g2] and pow2[g2] complete the term. C++ accumulates with a conditional subtraction, correct because total and term are already reduced; Python and Java take a remainder.

Only the C++ file carries checks. brute_count(n) applies the rules literally: a depth-first search over the bit mask of occupied seats, memoised in an unordered_map keyed by that mask, which finds the smallest neighbour count best, collects every empty seat attaining it and sums the counts of the resulting states. run_checkpoints then verifies $T(4)=8$ by that search and $T(10)=61632$ and the value of $T(1000)$ modulo $P$ by the formula, and main runs these checks unless it is given --skip-checkpoints. ^[9] The Python and Java ports contain the summation only, with no brute force and no checkpoints; Python also returns its total as a decimal string rather than as a number. In all three, the factorial and inverse-factorial tables run up to $n+1$ and the table of powers of two up to $n$. The surviving terms at $n=10^6$ only need factorial indices up to $500\,000$ and power indices up to $333\,333$, but they cannot simply be shortened: the inverse table is seeded from fact[n+1], and the loop loads fact[m] and fact[m-1] for every $m$ up to $n$, before the range test rejects a pair.

Complexity and Verification

Building the three tables costs $O(N)$ time and memory, the single modular exponentiation $O(\log P)$ multiplications, and the main loop $3N$ iterations of constant cost, so the method is $O(N)$ in time and memory; at $N=10^6$ exactly $500\,000$ of those iterations produce a term. The C++ and Java tables hold $3N+5$ signed 64-bit entries, $24\,000\,040$ bytes of elements at $N=10^6$, while a Python list stores references to arbitrary-precision integers instead of packed machine words and needs several times as much. On the machine used here, with Apple clang at -O2, CPython 3.11 and OpenJDK 17, C++ took $0.05$ seconds and peaked at $25$ MB, Java $0.05$ seconds and $70$ MB including its runtime, and Python $0.46$ seconds and $134$ MB; such figures depend on the hardware. The literal search is no route to the target, since the number of reachable masks of occupied seats grows exponentially; it serves as a small-case check, here over $1\le N\le20$.

These facts were checked by computation for this article. The literal search and the formula agree exactly, without any modulus, for every $N$ from $1$ to $20$: $1, 2, 4, 8, 32$, then $96, 456, 2016, 11232, 61632$, up to $T(20)=143113833676800$. These are the terms of OEIS A192008, the same rule phrased with phone booths, whose published formula becomes the sum of Section 9 when its $k$ is our $g_2$ and its own $m$ our $m-1-g_2$, the number of short internal gaps. ^[4] The payphone sequence A095236, where each newcomer takes one of the most distant free phones instead, runs $1, 2, 4, 8, 16, 36$ and already differs at $N=5$, where $T(5)=32$. ^[5]

Three further checks confirm the separate parts of the argument. Across the sizes $1\le N\le14$ there are $196$ phase-1 sets in all, and enumerating them showed that each has $m!$ phase-1 orders and $(s+g_2)!\,2^{g_2}(m-1)!$ completions, as Sections 6 to 8 claim. For every $N\le16$ the skeleton counts $w_s\binom{m-1}{g_2}$ add up to the number of maximal independent sets of the path found by direct enumeration; after the ten values of Section 5 these continue $21, 28, 37, 49, 65, 86$, with consecutive ratios approaching the $1.32471795\ldots$ of Section 5. ^[3] Finally the C++ checkpoints pass, and the three programs produce the same output at $N=10^6$, as does an independent implementation that sums the terms in a different order and replaces the powers-of-two table by a modular exponentiation; the outputs were compared byte by byte, never displayed.

Footnotes and References

The references below supply the background facts used above. The phase decomposition, the gap arithmetic, the counting of each phase and the operation counts are derived explicitly in this article.

Project Euler 364 — Comfortable Distance. The official statement puts N people one after another into a row of N seats under the three rules quoted above, defines T(N) as the number of possibilities, shows T(4) = 8 in a figure, states T(10) = 61632 and T(1000) mod 100000007 = 47255094, and asks for T(1000000) mod 100000007.
Wikipedia — Maximal independent set. It defines a maximal independent set as an independent set contained in no larger one, states that a vertex set is a maximal independent set if and only if it is an independent dominating set, which is the pair of properties proved in Section 3, and notes that the number of maximal independent sets of an n-vertex path graph is given by the Padovan sequence, with growth about 1.324718 per vertex.
OEIS A000931 — Padovan sequence. The entry defines a(n) = a(n−2) + a(n−3) with a(0) = 1 and a(1) = a(2) = 0, records that for n greater than 6 the term a(n) counts the maximal independent vertex sets of the (n−6)-path graph, and gives the asymptotics a(n) ~ r^n/(2r+3) with r = 1.3247179572447…, the real root of x³ = x + 1.
OEIS A192008 — modified linear phone booth sequence. It counts the ways to occupy n phone booths in a row one by one, each person picking a booth adjacent to the smallest number of occupied booths, which is the rule of this problem; its terms begin 1, 2, 4, 8, 32, 96, 456, 2016, 11232, 61632, and its formula sums (m+k+1)!·binomial(m+k, m)·2^k·(k+v1+v2)!·(m+k)! over v1, v2 in {0, 1} and nonnegative m, k with 2m + 3k = n − 1 − v1 − v2; its v1 + v2 is our s.
OEIS A095236 — payphone problem. Here each newcomer takes one of the phones most distant from those in use, the middle of a largest free span, with a span at the end of the row counted as having length 2L − 1; the terms 1, 2, 4, 8, 16, 36, 136 show that this maximum-distance rule is not the minimum-neighbour rule of this problem.
Wikipedia — Fermat's little theorem. For a prime p and any integer a the difference a^p − a is a multiple of p, and if a is coprime to p this is equivalent to a^(p−1) ≡ 1 (mod p). Section 11 draws the inverse a^(p−2) from that congruence; the article itself does not state that corollary.
Wikipedia — Exponentiation by squaring. The method removes the least significant binary digit of the exponent at each step, so it needs about log₂ n squarings plus one multiplication per one-bit, and the article names modular arithmetic as its typical use. This is exactly what mod_pow does.
cp-algorithms — Binomial coefficients. The page precomputes all factorials modulo m, notes that the binomial coefficient costs constant time once the inverse factorials are known, and gives the recurrence (x!)^(−1) ≡ ((x−1)!)^(−1)·x^(−1); rearranged as ((x−1)!)^(−1) ≡ (x!)^(−1)·x it is the downward pass that lets one modular inverse produce the whole table, as the code does.
C++ — Euler364.cpp. The linked immutable revision contains mod_pow, comfortable_arrangements_mod, brute_count and run_checkpoints, and is the implementation record for this article.

Problem 364 source code

C++

#include <cstdint>
#include <functional>
#include <iostream>
#include <string>
#include <unordered_map>
#include <vector>

namespace {

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

constexpr i64 kMod = 100000007LL;

i64 mod_pow(i64 base, i64 exp) {
    i64 result = 1 % kMod;
    base %= kMod;
    while (exp > 0) {
        if ((exp & 1LL) != 0LL) {
            result = static_cast<i64>((static_cast<__int128>(result) * base) % kMod);
        }
        base = static_cast<i64>((static_cast<__int128>(base) * base) % kMod);
        exp >>= 1LL;
    }
    return result;
}

i64 comfortable_arrangements_mod(const int n) {
    std::vector<i64> fact(static_cast<std::size_t>(n + 2), 1);
    for (int i = 1; i <= n + 1; ++i) {
        fact[static_cast<std::size_t>(i)] =
            static_cast<i64>((static_cast<__int128>(fact[static_cast<std::size_t>(i - 1)]) * i) % kMod);
    }

    std::vector<i64> inv_fact(static_cast<std::size_t>(n + 2), 1);
    inv_fact[static_cast<std::size_t>(n + 1)] = mod_pow(fact[static_cast<std::size_t>(n + 1)], kMod - 2);
    for (int i = n + 1; i >= 1; --i) {
        inv_fact[static_cast<std::size_t>(i - 1)] =
            static_cast<i64>((static_cast<__int128>(inv_fact[static_cast<std::size_t>(i)]) * i) % kMod);
    }

    std::vector<i64> pow2(static_cast<std::size_t>(n + 1), 1);
    for (int i = 1; i <= n; ++i) {
        pow2[static_cast<std::size_t>(i)] = (pow2[static_cast<std::size_t>(i - 1)] * 2LL) % kMod;
    }

    const int boundary_sum[3] = {0, 1, 2};
    const int boundary_mult[3] = {1, 2, 1};

    i64 total = 0;
    for (int m = 1; m <= n; ++m) {
        const i64 fm = fact[static_cast<std::size_t>(m)];
        const i64 fm1 = fact[static_cast<std::size_t>(m - 1)];

        for (int t = 0; t < 3; ++t) {
            const int s = boundary_sum[t];
            const int mult = boundary_mult[t];

            const int g2 = n - (2 * m - 1 + s);
            if (g2 < 0 || g2 > (m - 1)) {
                continue;
            }

            i64 comb = fm1;
            comb = static_cast<i64>((static_cast<__int128>(comb) * inv_fact[static_cast<std::size_t>(g2)]) % kMod);
            comb = static_cast<i64>(
                (static_cast<__int128>(comb) * inv_fact[static_cast<std::size_t>(m - 1 - g2)]) % kMod);

            i64 term = (static_cast<i64>(mult) * comb) % kMod;
            term = static_cast<i64>((static_cast<__int128>(term) * fm) % kMod);
            term = static_cast<i64>((static_cast<__int128>(term) * fm1) % kMod);
            term = static_cast<i64>((static_cast<__int128>(term) * fact[static_cast<std::size_t>(s + g2)]) % kMod);
            term = static_cast<i64>((static_cast<__int128>(term) * pow2[static_cast<std::size_t>(g2)]) % kMod);

            total += term;
            if (total >= kMod) {
                total -= kMod;
            }
        }
    }

    return total;
}

u64 brute_count(const int n) {
    const u64 full_mask = (n == 64) ? ~0ULL : ((1ULL << n) - 1ULL);
    std::unordered_map<u64, u64> memo;
    memo.reserve(1U << n);

    std::function<u64(u64)> dfs = [&](const u64 mask) -> u64 {
        if (mask == full_mask) {
            return 1ULL;
        }
        const auto it = memo.find(mask);
        if (it != memo.end()) {
            return it->second;
        }

        int best = 3;
        std::vector<int> choices;
        for (int i = 0; i < n; ++i) {
            if ((mask & (1ULL << i)) != 0ULL) {
                continue;
            }
            int occ = 0;
            if (i > 0 && (mask & (1ULL << (i - 1))) != 0ULL) {
                ++occ;
            }
            if (i + 1 < n && (mask & (1ULL << (i + 1))) != 0ULL) {
                ++occ;
            }
            if (occ < best) {
                best = occ;
                choices.clear();
                choices.push_back(i);
            } else if (occ == best) {
                choices.push_back(i);
            }
        }

        u64 total = 0;
        for (const int seat : choices) {
            total += dfs(mask | (1ULL << seat));
        }
        memo.emplace(mask, total);
        return total;
    };

    return dfs(0ULL);
}

bool run_checkpoints() {
    if (brute_count(4) != 8ULL) {
        std::cerr << "Checkpoint failed: T(4)\n";
        return false;
    }

    if (comfortable_arrangements_mod(10) != 61632LL) {
        std::cerr << "Checkpoint failed: T(10)\n";
        return false;
    }

    if (comfortable_arrangements_mod(1000) != 47255094LL) {
        std::cerr << "Checkpoint failed: T(1000) mod 100000007\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << comfortable_arrangements_mod(1000000) << '\n';
    return 0;
}

Python

def mod_pow(base, exp, mod):
    return pow(base, exp, mod)

def comfortable_arrangements_mod(n):
    mod = 100000007
    fact = [1] * (n + 2)
    for i in range(1, n + 2):
        fact[i] = (fact[i - 1] * i) % mod
        
    inv_fact = [1] * (n + 2)
    inv_fact[n + 1] = mod_pow(fact[n + 1], mod - 2, mod)
    for i in range(n + 1, 0, -1):
        inv_fact[i - 1] = (inv_fact[i] * i) % mod
        
    pow2 = [1] * (n + 1)
    for i in range(1, n + 1):
        pow2[i] = (pow2[i - 1] * 2) % mod
        
    boundary_sum = [0, 1, 2]
    boundary_mult = [1, 2, 1]
    
    total = 0
    for m in range(1, n + 1):
        fm = fact[m]
        fm1 = fact[m - 1]
        
        for t in range(3):
            s = boundary_sum[t]
            mult = boundary_mult[t]
            
            g2 = n - (2 * m - 1 + s)
            if g2 < 0 or g2 > (m - 1):
                continue
                
            comb = fm1
            comb = (comb * inv_fact[g2]) % mod
            comb = (comb * inv_fact[m - 1 - g2]) % mod
            
            term = (mult * comb) % mod
            term = (term * fm) % mod
            term = (term * fm1) % mod
            term = (term * fact[s + g2]) % mod
            term = (term * pow2[g2]) % mod
            
            total = (total + term) % mod
            
    return str(total)

def solve():
    return comfortable_arrangements_mod(1000000)

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

Java

public class Euler364 {
    static long modPow(long base, long exp, long mod) {
        long result = 1 % mod;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                result = (result * base) % mod;
            base = (base * base) % mod;
            exp >>= 1;
        }
        return result;
    }

    static long comfortableArrangementsMod(int n) {
        long mod = 100000007L;
        long[] fact = new long[n + 2];
        fact[0] = 1;
        for (int i = 1; i <= n + 1; ++i) {
            fact[i] = (fact[i - 1] * i) % mod;
        }

        long[] invFact = new long[n + 2];
        invFact[n + 1] = modPow(fact[n + 1], mod - 2, mod);
        for (int i = n + 1; i >= 1; --i) {
            invFact[i - 1] = (invFact[i] * i) % mod;
        }

        long[] pow2 = new long[n + 1];
        pow2[0] = 1;
        for (int i = 1; i <= n; ++i) {
            pow2[i] = (pow2[i - 1] * 2) % mod;
        }

        int[] boundarySum = { 0, 1, 2 };
        int[] boundaryMult = { 1, 2, 1 };

        long total = 0;
        for (int m = 1; m <= n; ++m) {
            long fm = fact[m];
            long fm1 = fact[m - 1];

            for (int t = 0; t < 3; ++t) {
                int s = boundarySum[t];
                int mult = boundaryMult[t];

                int g2 = n - (2 * m - 1 + s);
                if (g2 < 0 || g2 > (m - 1))
                    continue;

                long comb = fm1;
                comb = (comb * invFact[g2]) % mod;
                comb = (comb * invFact[m - 1 - g2]) % mod;

                long term = (mult * comb) % mod;
                term = (term * fm) % mod;
                term = (term * fm1) % mod;
                term = (term * fact[s + g2]) % mod;
                term = (term * pow2[g2]) % mod;

                total = (total + term) % mod;
            }
        }

        return total;
    }

    public static String solve() {
        return String.valueOf(comfortableArrangementsMod(1000000));
    }

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