Problem 837: Amidakuji

Project Euler Problem 837 Solution

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

Problem Summary There are three vertical lines and two allowed adjacent swaps. Operation $A$ swaps the first and second lines, and operation $B$ swaps the second and third lines. Among all words that contain exactly $m$ copies of $A$ and exactly $n$ copies of $B$, we must count those whose total permutation is the identity. The required answer is taken modulo $M=1234567891$. A direct dynamic program over all prefixes and all six permutations is useful only for tiny checks. For the real input sizes, the key is to compress the word into consecutive pairs and count those pairs analytically. Mathematical Approach Step 1: Model the Ladder as Permutations Work in the symmetric group $S_3$: $$A=(12),\qquad B=(23).$$ Both generators are transpositions, so they are odd permutations. The identity permutation is even....

Detailed mathematical approach

Problem Summary

There are three vertical lines and two allowed adjacent swaps. Operation $A$ swaps the first and second lines, and operation $B$ swaps the second and third lines. Among all words that contain exactly $m$ copies of $A$ and exactly $n$ copies of $B$, we must count those whose total permutation is the identity.

The required answer is taken modulo $M=1234567891$. A direct dynamic program over all prefixes and all six permutations is useful only for tiny checks. For the real input sizes, the key is to compress the word into consecutive pairs and count those pairs analytically.

Mathematical Approach

Step 1: Model the Ladder as Permutations

Work in the symmetric group $S_3$:

$$A=(12),\qquad B=(23).$$

Both generators are transpositions, so they are odd permutations. The identity permutation is even. Therefore, if $m+n$ is odd, then no word can evaluate to the identity and

$$F(m,n)=0.$$

From now on, assume $m+n$ is even and write

$$m+n=2t,\qquad t=\frac{m+n}{2}.$$

Step 2: Group the Word into Consecutive Pairs

Split any word of length $2t$ into $t$ adjacent blocks:

$$\bigl(x_1x_2\bigr)\bigl(x_3x_4\bigr)\cdots\bigl(x_{2t-1}x_{2t}\bigr),\qquad x_i\in\{A,B\}.$$

Now set

$$e=\text{identity},\qquad r=AB=(123),\qquad r^2=BA=(132).$$

Each pair belongs to the cyclic subgroup $\{e,r,r^2\}$, because

$$AA=e,\qquad BB=e,\qquad AB=r,\qquad BA=r^2.$$

Suppose exactly $k$ of the $t$ pairs are mixed, meaning they are either $AB$ or $BA$. Then the remaining pairs must be neutral pairs $AA$ and $BB$. Their counts are forced:

$$u=\frac{m-k}{2},\qquad v=\frac{n-k}{2}.$$

So $k$ is admissible exactly when

$$0\le k\le \min(m,n),\qquad k\equiv m\pmod 2.$$

Once $k$ is fixed, the number of ways to choose which of the $t$ pair positions are $AA$, which are $BB$, and which are mixed is the multinomial coefficient

$$T_k=\frac{t!}{u!\,v!\,k!}=\frac{t!}{\left(\frac{m-k}{2}\right)!\left(\frac{n-k}{2}\right)!k!}.$$

This counts the block layout, but it does not yet distinguish whether a mixed block is $AB$ or $BA$.

Step 3: Count the Mixed Blocks That Multiply to the Identity

Inside the $k$ mixed positions, let exactly $j$ blocks be $AB=r$ and the remaining $k-j$ blocks be $BA=r^2$. Their product is

$$r^j(r^2)^{k-j}=r^{j+2(k-j)}=r^{2j-k}.$$

This equals the identity exactly when

$$2j-k\equiv 0\pmod 3.$$

Therefore the number of valid orientations of the $k$ mixed blocks is

$$R_k=\sum_{\substack{0\le j\le k\\2j-k\equiv 0\pmod 3}}\binom{k}{j}.$$

This sum has a clean closed form. Let $\omega\neq 1$ be a complex cube root of unity, so $\omega^3=1$. A roots-of-unity filter gives

$$R_k=\frac{1}{3}\sum_{q=0}^{2}\omega^{-2kq}(1+\omega^q)^k.$$

Evaluating the three terms yields

$$R_k=\frac{2^k+2(-1)^k}{3}.$$

Step 4: Combine the Layout Count and the Orientation Count

For each admissible $k$, the choices of block layout and the choices of mixed-block orientations are independent. Hence

$$F(m,n)=\sum_{\substack{0\le k\le \min(m,n)\\k\equiv m\pmod 2}} T_k\,R_k \pmod{M}.$$

Substituting the explicit formulas gives

$$F(m,n)=\sum_{\substack{0\le k\le \min(m,n)\\k\equiv m\pmod 2}}\frac{t!}{\left(\frac{m-k}{2}\right)!\left(\frac{n-k}{2}\right)!k!}\cdot\frac{2^k+2(-1)^k}{3}\pmod{M}.$$

This is the closed summation used by the implementations.

Step 5: Derive the Recurrence Used by the Implementation

Let $k_0=m\bmod 2$, the smallest admissible value of $k$. The initial combinatorial factor is

$$T_{k_0}= \begin{cases} \binom{t}{m/2}, & k_0=0,\\ t\binom{t-1}{(m-1)/2}, & k_0=1. \end{cases}$$

For later terms, taking the ratio of consecutive admissible values gives

$$\frac{T_{k+2}}{T_k}=\frac{(m-k)(n-k)}{4(k+1)(k+2)}.$$

So once one term is known, the next one is obtained with constant-time modular arithmetic. Also, when $k$ increases by $2$, the factor $2^k$ is multiplied by $4$, while $(-1)^k$ stays unchanged because the parity of $k$ never changes inside the loop.

Worked Example: $m=n=3$

Here $m+n=6$, so $t=3$, and admissible values of $k$ are $1$ and $3$.

For $k=1$, we have

$$T_1=\frac{3!}{1!\,1!\,1!}=6,\qquad R_1=\frac{2^1+2(-1)^1}{3}=0.$$

So all layouts with exactly one mixed pair contribute nothing, because a single nontrivial element of order $3$ cannot be the identity.

For $k=3$, we get

$$T_3=\frac{3!}{0!\,0!\,3!}=1,\qquad R_3=\frac{2^3+2(-1)^3}{3}=2.$$

Thus

$$F(3,3)=1\cdot 2=2.$$

The two words are the fully alternating ones:

$$ABABAB,\qquad BABABA.$$

How the Code Works

The C++, Python, and Java implementations all evaluate the same formula modulo $1234567891$. They first handle the parity test $m+n$ odd $\Rightarrow 0$, then set $t=(m+n)/2$ and start from the smallest admissible $k$.

The initial factor $T_{k_0}$ is computed through a binomial coefficient modulo $M$. Since modular division is needed, the implementation uses modular inverses under the prime modulus, obtained by fast exponentiation. To avoid computing one inverse at a time, it builds inverses for long consecutive ranges in batches and reuses them inside the product.

After that, the program walks through $k,k+2,k+4,\dots$. At each step it updates the current multinomial factor with the ratio

$$\frac{(m-k)(n-k)}{4(k+1)(k+2)}$$

updates the power $2^k$ by multiplying by $4$, multiplies by the modular inverse of $3$, and adds the current contribution to the running answer. The whole computation therefore uses the closed form directly and never builds a large $O(mn)$ dynamic-programming table over the six states of $S_3$.

Complexity Analysis

The number of admissible values of $k$ is $1+\left\lfloor\frac{\min(m,n)-k_0}{2}\right\rfloor$, so the main summation is linear in $\min(m,n)$. Computing the starting binomial factor is also linear in the smaller half-count, so the total running time is $O(\min(m,n))$ modular multiplications.

The implementations store only a few scalar values plus a temporary batch of modular inverses. If the batch size is denoted by $B$, the extra memory is $O(B)$; with the batch size treated as a fixed engineering constant, the working memory is effectively constant with respect to $m$ and $n$.

Footnotes and References

Problem page: Project Euler 837
Symmetric group: Wikipedia — Symmetric group
Roots of unity: Wikipedia — Root of unity
Multinomial theorem: Wikipedia — Multinomial theorem
Fermat's little theorem: Wikipedia — Fermat's little theorem

Problem 837 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

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

static constexpr u64 kMod = 1'234'567'891ULL;
static constexpr u64 kInv3 = 823'045'261ULL;
static constexpr u64 kBlock = 8'000'000ULL;
static constexpr u64 kBarrett = static_cast<u64>((static_cast<u128>(1) << 64) / kMod);

static inline u64 mod_reduce(u64 x) {
    const u64 q = static_cast<u64>((static_cast<u128>(x) * kBarrett) >> 64);
    u64 r = x - q * kMod;
    if (r >= kMod) r -= kMod;
    if (r >= kMod) r -= kMod;
    return r;
}

static inline u64 mod_mul(u64 a, u64 b) {
    return mod_reduce(a * b);
}

static u64 mod_pow(u64 base, u64 exp) {
    u64 out = 1;
    base %= kMod;
    while (exp > 0) {
        if (exp & 1ULL) out = mod_mul(out, base);
        base = mod_mul(base, base);
        exp >>= 1ULL;
    }
    return out;
}

static inline u64 mod_inv(u64 x) {
    return mod_pow(x, kMod - 2);
}

static void inverse_range(u64 start, u64 len, std::vector<u64>& pref) {
    pref.resize(static_cast<std::size_t>(len));
    u64 cur = start % kMod;
    pref[0] = cur;
    for (u64 i = 1; i < len; ++i) {
        if (++cur == kMod) cur = 0;
        pref[static_cast<std::size_t>(i)] = mod_mul(pref[static_cast<std::size_t>(i - 1)], cur);
    }

    u64 inv_prod = mod_inv(pref[static_cast<std::size_t>(len - 1)]);
    cur = (start + len - 1) % kMod;
    for (u64 i = len; i-- > 0;) {
        const u64 prev = (i == 0) ? 1 : pref[static_cast<std::size_t>(i - 1)];
        pref[static_cast<std::size_t>(i)] = mod_mul(inv_prod, prev);
        inv_prod = mod_mul(inv_prod, cur);
        if (cur == 0) cur = kMod - 1;
        else --cur;
    }
}

static u64 binom_mod(u64 n, u64 k) {
    if (k > n) return 0;
    k = std::min(k, n - k);
    if (k == 0) return 1;

    u64 out = 1;
    const u64 offset = n - k;
    std::vector<u64> invs;
    invs.reserve(static_cast<std::size_t>(kBlock));
    for (u64 l = 1; l <= k; l += kBlock) {
        const u64 r = std::min(k, l + kBlock - 1);
        const u64 len = r - l + 1;
        inverse_range(l, len, invs);
        u64 num = (offset + l) % kMod;
        for (u64 i = 0; i < len; ++i) {
            out = mod_mul(out, num);
            out = mod_mul(out, invs[static_cast<std::size_t>(i)]);
            if (++num == kMod) num = 0;
        }
    }
    return out;
}

static u64 solve(u64 m, u64 n) {
    if ((m + n) & 1ULL) return 0;

    const u64 t = (m + n) >> 1;
    const u64 parity = m & 1ULL;
    const u64 max_k = std::min(m, n);

    u64 k = parity;
    u64 a = (m - parity) >> 1;
    u64 b = (n - parity) >> 1;

    u64 weight = 0;
    if (parity == 0ULL) {
        weight = binom_mod(t, a);
    } else {
        weight = mod_mul(t % kMod, binom_mod(t - 1ULL, a));
    }

    u64 pow2k = (parity == 0ULL) ? 1ULL : 2ULL;
    const u64 sign_part = (parity == 0ULL) ? 2ULL : (kMod - 2ULL);  // 2*(-1)^k

    u64 answer = 0ULL;
    auto add_current = [&]() {
        u64 f3 = pow2k + sign_part;
        if (f3 >= kMod) f3 -= kMod;
        f3 = mod_mul(f3, kInv3);
        const u64 inc = mod_mul(weight, f3);
        answer += inc;
        if (answer >= kMod) answer -= kMod;
    };

    add_current();
    if (k == max_k) return answer;

    u64 updates = (max_k - k) >> 1;
    std::vector<u64> invs;
    invs.reserve(static_cast<std::size_t>(2 * kBlock));
    u64 a_mod = a % kMod;
    u64 b_mod = b % kMod;

    while (updates > 0ULL) {
        const u64 len = std::min<u64>(updates, kBlock);
        inverse_range(k + 1ULL, 2ULL * len, invs);

        for (u64 i = 0; i < len; ++i) {
            const u64 inv_k1 = invs[static_cast<std::size_t>(2ULL * i)];
            const u64 inv_k2 = invs[static_cast<std::size_t>(2ULL * i + 1ULL)];

            weight = mod_mul(weight, a_mod);
            weight = mod_mul(weight, b_mod);
            weight = mod_mul(weight, inv_k1);
            weight = mod_mul(weight, inv_k2);

            --a;
            --b;
            if (a_mod == 0) a_mod = kMod - 1;
            else --a_mod;
            if (b_mod == 0) b_mod = kMod - 1;
            else --b_mod;
            k += 2ULL;
            pow2k = mod_mul(pow2k, 4ULL);

            add_current();
        }
        updates -= len;
    }

    return answer;
}

static u64 brute_small(int m, int n) {
    const std::array<std::array<int, 3>, 6> perms = {{
        {{0, 1, 2}}, {{0, 2, 1}}, {{1, 0, 2}},
        {{1, 2, 0}}, {{2, 0, 1}}, {{2, 1, 0}}
    }};

    std::array<int, 6> trans_ab{};
    std::array<int, 6> trans_bc{};
    for (int i = 0; i < 6; ++i) {
        auto p1 = perms[i];
        std::swap(p1[0], p1[1]);
        auto p2 = perms[i];
        std::swap(p2[1], p2[2]);
        for (int j = 0; j < 6; ++j) {
            if (perms[j] == p1) trans_ab[i] = j;
            if (perms[j] == p2) trans_bc[i] = j;
        }
    }

    std::vector<std::vector<std::array<u64, 6>>> dp(
        m + 1, std::vector<std::array<u64, 6>>(n + 1, {0, 0, 0, 0, 0, 0}));
    dp[0][0][0] = 1;

    for (int i = 0; i <= m; ++i) {
        for (int j = 0; j <= n; ++j) {
            if (i == 0 && j == 0) continue;
            std::array<u64, 6> cur = {0, 0, 0, 0, 0, 0};
            if (i > 0) {
                for (int s = 0; s < 6; ++s) {
                    cur[trans_ab[s]] += dp[i - 1][j][s];
                }
            }
            if (j > 0) {
                for (int s = 0; s < 6; ++s) {
                    cur[trans_bc[s]] += dp[i][j - 1][s];
                }
            }
            dp[i][j] = cur;
        }
    }
    return dp[m][n][0];
}

int main() {
    for (int m = 0; m <= 8; ++m) {
        for (int n = 0; n <= 8; ++n) {
            assert(solve(m, n) == brute_small(m, n) % kMod);
        }
    }
    assert(solve(3, 3) == 2);
    assert(solve(123, 321) == 172'633'303ULL);

    std::cout << solve(123'456'789ULL, 987'654'321ULL) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

public class Euler837 {
    static final long kMod = 1234567891L;
    static final long kInv3 = 823045261L;
    static final int kBlock = 8000000;

    static long modPow(long base, long exp) {
        long out = 1;
        base %= kMod;
        while (exp > 0) {
            if ((exp & 1) == 1)
                out = (out * base) % kMod;
            base = (base * base) % kMod;
            exp >>= 1;
        }
        return out;
    }

    static long modInv(long x) {
        return modPow(x, kMod - 2);
    }

    static long[] inverseRange(long start, int len) {
        if (len == 0)
            return new long[0];
        long[] pref = new long[len];
        long cur = start % kMod;
        pref[0] = cur;
        for (int i = 1; i < len; ++i) {
            if (++cur == kMod)
                cur = 0;
            pref[i] = (pref[i - 1] * cur) % kMod;
        }

        long invProd = modInv(pref[len - 1]);
        cur = (start + len - 1) % kMod;

        for (int i = len - 1; i >= 0; --i) {
            long prev = (i == 0) ? 1 : pref[i - 1];
            pref[i] = (invProd * prev) % kMod;
            invProd = (invProd * cur) % kMod;
            if (cur == 0)
                cur = kMod - 1;
            else
                --cur;
        }
        return pref;
    }

    static long binomMod(long n, long k) {
        if (k > n)
            return 0;
        k = Math.min(k, n - k);
        if (k == 0)
            return 1;

        long out = 1;
        long offset = n - k;

        for (long l = 1; l <= k; l += kBlock) {
            long r = Math.min(k, l + kBlock - 1);
            int len = (int) (r - l + 1);
            long[] invs = inverseRange(l, len);
            long num = (offset + l) % kMod;
            for (int i = 0; i < len; ++i) {
                out = (out * num) % kMod;
                out = (out * invs[i]) % kMod;
                if (++num == kMod)
                    num = 0;
            }
        }
        return out;
    }

    static long solve(long m, long n) {
        if (((m + n) & 1) == 1)
            return 0;

        long t = (m + n) >> 1;
        long parity = m & 1;
        long maxK = Math.min(m, n);

        long k = parity;
        long a = (m - parity) >> 1;
        long b = (n - parity) >> 1;

        long weight = 0;
        if (parity == 0) {
            weight = binomMod(t, a);
        } else {
            weight = ((t % kMod) * binomMod(t - 1, a)) % kMod;
        }

        long pow2k = (parity == 0) ? 1 : 2;
        long signPart = (parity == 0) ? 2 : (kMod - 2);

        long answer = 0;

        long f3 = pow2k + signPart;
        if (f3 >= kMod)
            f3 -= kMod;
        f3 = (f3 * kInv3) % kMod;
        long inc = (weight * f3) % kMod;
        answer = (answer + inc) % kMod;

        if (k == maxK)
            return answer;

        long updates = (maxK - k) >> 1;
        long aMod = a % kMod;
        long bMod = b % kMod;

        while (updates > 0) {
            int len = (int) Math.min(updates, kBlock);
            long[] invs = inverseRange(k + 1, 2 * len);

            for (int i = 0; i < len; ++i) {
                long invK1 = invs[2 * i];
                long invK2 = invs[2 * i + 1];

                weight = (weight * aMod) % kMod;
                weight = (weight * bMod) % kMod;
                weight = (weight * invK1) % kMod;
                weight = (weight * invK2) % kMod;

                --a;
                --b;
                if (aMod == 0)
                    aMod = kMod - 1;
                else
                    --aMod;
                if (bMod == 0)
                    bMod = kMod - 1;
                else
                    --bMod;
                k += 2;
                pow2k = (pow2k * 4) % kMod;

                f3 = pow2k + signPart;
                if (f3 >= kMod)
                    f3 -= kMod;
                f3 = (f3 * kInv3) % kMod;
                inc = (weight * f3) % kMod;
                answer = (answer + inc) % kMod;
            }
            updates -= len;
        }

        return answer;
    }

    public static String solveStr() {
        return Long.toString(solve(123456789L, 987654321L));
    }

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