Problem 227: The Chase

Project Euler Problem 227 Solution

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

Problem Summary One hundred players sit evenly around a circle, and two dice start in opposite positions. In each round, each die independently moves one seat in one direction with probability $1/6$, stays where it is with probability $2/3$, or moves one seat in the other direction with probability $1/6$. The chase ends as soon as both dice are in the same player's hands. Because the initial separation is 50 seats on a 100-player circle, the required quantity is the expected hitting time of the meeting state starting from circular distance $50$. The implementations evaluate this expectation as $3780.6186217845$ rounds. Mathematical Approach The key simplification is to ignore absolute seat numbers and keep only the circular distance between the two dice. That turns the game into a finite absorbing Markov chain whose transition law can be written down exactly. Rotational invariance collapses the state space Let $k\in\{0,1,\dots,n-1\}$ be the clockwise distance from the first die to the second on a circle of $n=100$ players. If both dice are rotated together by the same number of seats, nothing about future play changes, so the expectation depends only on $k$, not on the absolute positions. Define $E_k$ as the expected number of remaining rounds until the dice meet when the current distance is $k$....

Detailed mathematical approach

Problem Summary

An even number $n$ of players sit around a table, and two opposite players start with one die each. On every turn, the two players holding a die roll it: a $1$ passes the die to the neighbour on the left, a $6$ to the neighbour on the right, and any other number lets the player keep it. The game ends when, after the rolls and passes of a turn, one player holds both dice; that player loses. For $n=100$ we need the expected number of turns, rounded to ten significant digits; no values for smaller tables are given. ^[1]

The whole solution in one sentence: only the number of seats from one die to the other matters, this distance performs a random walk around the table that stops at $0$, and the expected remaining lengths from the $99$ nonzero distances are the unique solution of $99$ linear equations, which the program solves by elimination.

Mathematical Approach

We reduce the dice to one number, set up equations for the expected length, prove that they determine it, solve small cases by hand, and check how accurately a computer solves them.

1. One turn of one die

Number the seats $0,1,\dots,n-1$ so that a pass to the left adds $1$ to the seat number of the die and a pass to the right subtracts $1$, modulo $n$; the other choice only mirrors the game. In one turn the seat of a die changes by $X=+1$ or $X=-1$ with probability $1/6$ each and by $X=0$ with probability $2/3$, independently of the other die and of earlier turns.

2. Only the distance matters

Call the dice $A$ and $B$, and let $D$ be the seat number of $B$ minus that of $A$, modulo $n$. Turning the table changes nothing, so only $D$ matters. The game starts at $D=n/2$ and ends exactly when $D=0$. In a turn, $D$ changes by $Z=X_B-X_A$. Among the $36$ equally likely pairs of rolls:

$Z=+2$ in $1$ pair: $B$ shows $1$ and $A$ shows $6$;
$Z=+1$ in $8$ pairs: $B$ shows $1$ and $A$ one of $2,\dots,5$, or $B$ shows one of $2,\dots,5$ and $A$ shows $6$;
$Z=0$ in $18$ pairs: both dice stay ($16$ pairs) or both move the same way ($2$ pairs);
$Z=-1$ in $8$ pairs and $Z=-2$ in $1$ pair, the mirror images of the first two cases.

So the probability $p_z$ of the change $z$ is

$$p_{\pm2}=\frac1{36},\quad p_{\pm1}=\frac29,\quad p_0=\frac12.$$

The change does not depend on the past, so $D$ is a Markov chain: a random walk on the cycle of residues modulo $n$ that is absorbed at $0$. ^[2]

3. The dice can pass each other

A change of $2$ lets the dice swap places without meeting. With $A$ on seat $0$ and $B$ on seat $1$, the rolls $1$ for $A$ and $6$ for $B$ put $A$ on seat $1$ and $B$ on seat $0$: nobody holds both dice, and the new distance is $1-2\equiv n-1$. From $D=2$ the same rolls bring both dice to the seat between them and end the game. Adding $Z$ modulo $n$ and stopping only at $D=0$ handles both cases. For $n\ge6$, a turn ends the game with probability $2/9$ from $D\equiv\pm1$, with probability $1/36$ from $D\equiv\pm2$, and never from the other distances.

4. Equations for the expected length

Start with two players, where the only nonzero distance is $1$. The changes $\pm1$ end the game, while $-2$, $0$ and $+2$ lead back to distance $1$, together with probability $\frac1{36}+\frac12+\frac1{36}=\frac59$. The first turn always counts, and with probability $5/9$ the game is then where it started, so the expected length $E_1$ satisfies

$$E_1=1+\frac59E_1,\qquad E_1=\frac94.$$

Indeed, exactly one die moves with probability $2\cdot\frac13\cdot\frac23=\frac49$ per turn, and $9/4$ is the reciprocal of that probability.

In general, let $E_d$ be the expected number of remaining turns at distance $d$; then $E_0=0$, and the answer is $E_{n/2}$. From $d\ne0$ one turn is played, the distance becomes $(d+z)\bmod n$ with probability $p_z$, and the game continues like a new game from there. Averaging over the five changes, which is called first-step analysis, gives for $d=1,\dots,n-1$

$$E_d=1+\sum_{z=-2}^{2}p_z\,E_{(d+z)\bmod n}.$$

With $Q$ the $(n-1)\times(n-1)$ matrix of one-turn probabilities between nonzero distances and $\mathbf 1$ the vector of ones, this reads $(I-Q)E=\mathbf 1$, the classical equation for absorbing Markov chains: the expected times to absorption are the row sums $N\mathbf 1$ of the fundamental matrix $N=(I-Q)^{-1}$. ^[2] The derivation assumed that every $E_d$ is finite; Section 5 proves it.

5. Every game ends

From any $d\ne0$, single steps the shorter way round, each of probability $2/9$, reach $0$ within $n/2$ turns. So each block of $n/2$ turns ends the game with probability at least $q=(2/9)^{n/2}$, whatever happened before, and the game lasts more than $kn/2$ turns with probability at most $(1-q)^k$. Hence the game ends with probability $1$, and

$$E_d\le\frac n2\sum_{k\ge0}(1-q)^k=\frac{n}{2q}.$$

The bound is enormous, but finiteness is all Section 4 needs. Likewise the row sums of $Q^k$, the probabilities that the game is still on after $k$ turns, are at most $(1-q)^{\lfloor 2k/n\rfloor}$, so $I+Q+Q^2+\cdots$ converges, and since $(I-Q)(I+\cdots+Q^k)=I-Q^{k+1}$, its sum is $N$.

6. The equations have exactly one solution

If $u$ and $v$ both solve the equations with $u_0=v_0=0$, then $w=u-v$ has $w_0=0$ and $w_d=\sum_z p_z\,w_{(d+z)\bmod n}$ for $d\ne0$: each $w_d$ is an average, with positive weights, of itself and its four neighbours. Suppose the largest entry $M$ of $w$ is positive, attained at some $d\ne0$. An average equals the largest of its values only if all of them are equal to it, so $w_{d-1}=M$; repeating this from $d-1$, $d-2$ and so on, we reach $w_0=M\gt0$, a contradiction. So $w\le0$, in the same way $w\ge0$, and $u=v$: the matrix $I-Q$ is invertible, and an exact solution of the equations gives the true expectations. Since $p_z=p_{-z}$, the vector with the entries $E_{n-d}$ solves the same equations, so uniqueness also gives the symmetry $E_d=E_{n-d}$.

7. Worked example: four players

For $n=4$, symmetry gives $E_3=E_1$. From $d=1$ the changes $-2,-1,0,+1,+2$ lead to the distances $3,0,1,2,3$, the first of them being the crossing of Section 3; from $d=2$ they lead to $0,1,2,3,0$. Hence

$$\begin{aligned}E_1&=1+\frac59E_1+\frac29E_2,\\E_2&=1+\frac49E_1+\frac12E_2,\end{aligned}$$

that is, $4E_1-2E_2=9$ and $9E_2-8E_1=18$. Doubling the first equation and adding the second gives $5E_2=36$, so $E_2=36/5=7.2$ and $E_1=(9+2E_2)/4=117/20$. The value $36/5$ and the $9/4$ of Section 4 are the checkpoints of the C++ program.

8. Larger tables

After using symmetry, $n=6$ leaves three unknowns and gives

$$\begin{aligned}E_1&=\frac{419}{44},\qquad E_2=\frac{152}{11},\\E_3&=\frac{675}{44}\approx15.34,\end{aligned}$$

and $n=8$ gives $E_4=1296/49\approx26.45$. The values grow roughly like $n^2$. For comparison, the drunkard's walk on $0,1,2,3,4$, which steps to either neighbour with probability $1/2$ and stops at $0$ and $4$, has the expected times $k(4-k)$, that is, $3$, $4$, $3$ from $k=1,2,3$. ^[2]

9. Solving the 99 equations

For $n=100$ the program keeps all $99$ unknowns, so the augmented matrix $[\,I-Q\mid\mathbf 1\,]$ has $99$ rows and $100$ columns. Gauss–Jordan elimination treats the columns in order. ^[3] It takes as pivot the entry of largest absolute value among the unused rows (partial pivoting), swaps that row into place, divides it by the pivot, and subtracts multiples of it from all other rows until the rest of the column is $0$. At the end the left part is the identity matrix and the last column holds $E_1,\dots,E_{99}$; a nonzero pivot always exists, because $I-Q$ is invertible.

The diagonal entry $1-p_0=1/2$ of each row is at least the sum of the absolute values of the other entries, strictly only for the distances $1$, $2$, $98$, $99$, and $Q$ is symmetric, so the columns are diagonally dominant too. Strict column dominance would make pivoting unnecessary; ^[4] ours is not strict everywhere, yet the pivot search never swaps two rows for any even $n\le100$ (we counted).

10. How many digits survive the rounding

Binary64 numbers, used by Python, Java and our C++ builds, keep about $16$ significant digits, and every operation rounds. The condition number $\kappa=\|A\|\,\|A^{-1}\|$ of $A=I-Q$ tells how many may be lost: roughly $k$ digits if $\kappa=10^k$. ^[5] In the maximum norm, the largest absolute row sum, $\|A\|_\infty=1$, since the rows away from distance $0$ have absolute row sum $1/2+1/2$ and the other four less. The inverse $N=I+Q+Q^2+\cdots$ is nonnegative with the row sums $N\mathbf 1=E$, so $\|A^{-1}\|_\infty=\max_d E_d$, which for $100$ players is $E_{50}$, as $E_d$ increases up to $d=50$ (checked on the exact solution). The condition number is thus the very expectation we compute. It is below $10^4$, so by the rule of thumb about four of the sixteen digits are at risk, besides the elimination's own rounding errors. That makes ten significant digits plausible; the exact comparison in the verification section confirms them.

How the Code Works

The C++ program ^[7] reads --players=N, default $100$, and --skip-checkpoints in parse_arguments; an unknown argument, an odd $N$ or $N\lt2$ gives exit code $1$. parse_int_after_prefix accumulates the digits in an int without overflow check, so an oversized count is undefined behaviour. Unless skipped, run_checkpoints requires solve(2) and solve(4) to be within $10^{-12}$ of $2.25$ and $7.2$ (Sections 4 and 7), or the exit code is $2$. main prints the result as a double with ten decimals; rounding to ten significant digits is left to the reader.

solve keeps the augmented matrix of Section 9 in a, $m=n-1$ rows of $m+1$ long double entries, and the changes and probabilities of Section 2 in delta and prob. Row k - 1 is distance $k$: its diagonal entry and right-hand side are set to $1$, and for each change the target nk ${}=(k+\delta)\bmod n$, made nonnegative because the C++ remainder can be negative, is skipped if it is $0$; otherwise the probability is subtracted in column nk - 1. With -=, changes with the same target add up: for $n=2$, the changes $-2$, $0$ and $+2$ all return to distance $1$, giving the diagonal entry $1-5/9=4/9$ of Section 4.

The elimination loop is the Gauss–Jordan elimination of Section 9: a pivot search with fabsl and a strict comparison, std::swap of two row vectors, division of the pivot row from column col on, and subtraction of factor times the pivot row from every other row. Rows whose factor is below $10^{-21}$ in absolute value are skipped; in the Python version this changes no bit of the result for $n=100$ (we compared). The answer is the last entry of row $n/2-1$.

Python's solve(players=100), called with its default, and Java's solve(), which fixes $100$ players, run the same loops on a list of float lists and on a double[][] array, without command-line arguments or checkpoints. Python's test nk < 0 never succeeds, as its % is never negative here; Java formats with Locale.US. Both then compare the ten-decimal string with one fixed string and, on a match, raise its last digit by one, which makes their output equal to that of the C++ program on our ARM Mac.

Python floats and Java doubles are IEEE binary64 numbers; ^[8] the C++ long double is 80-bit extended with GCC on x86 Linux, but plain double with Microsoft Visual C++ and on 64-bit ARM under Windows and macOS. ^[9] On our ARM Mac, Clang and GCC fuse the update a[row][j] -= factor * a[col][j] into one multiply–subtract that rounds once; ^[10] this alone raises the tenth printed decimal by one against unpatched Python and Java, and with -ffp-contract=off the C++ program prints their string. The 80-bit format named in the comment at the patch would, by our emulation, print the correctly rounded tenth decimal instead. Such differences of order $10^{-10}$ leave the ten significant digits unchanged, as the verification shows.

Complexity and Verification

Filling the $(n-1)n$ entries, $9900$ for $n=100$, takes $O(n^2)$ time and memory. For column $c=0,\dots,m-1$ the elimination updates $m+1-c$ entries in each of the other $m-1$ rows, at most $(m-1)\bigl((m+1)(m+2)/2-1\bigr)$ updates, $494\,802$ for $m=99$, so the time is $O(n^3)$; the skipped rows leave $181\,682$ updates, and each program runs in well under a tenth of a second on our machine. Keeping only the distances $1,\dots,n/2$, with each target $r$ replaced by $\min(r,n-r)$, gives a banded matrix, since a turn changes the shorter distance by at most $2$; a banded solver would need only $O(n)$ operations.

Simulation is no alternative. The second moments $S_d$ of the length satisfy $(I-Q)S=2E-\mathbf 1$, as squaring $T=1+T'$ and averaging as in Section 4 shows, and they give a standard deviation of about $0.82$ times the mean, for four players (variance $873/25$ at $d=2$) as for $100$. With the standard error $\sigma/\sqrt M$ of $M$ games, ^[6] a relative precision of $10^{-10}$ would take about $10^{20}$ games.

Besides the C++ checkpoints, we ran these checks, which are not part of the published programs:

exact solutions in Python fractions for every even $n\le40$ are symmetric, increase up to $d=n/2$ and agree in all ten printed decimals with the C++ output for --players=$n$;
$10^6$ simulated games of the real game, with two dice moving between seats and no distance model, for each of $n=2,4,6,8,10$ give averages within $1.1$ standard errors of the exact values, for example $7.1965\pm0.0059$ against $36/5$;
for $n=100$, the exact solution of the reduced system with $50$ unknowns matches the $99$ equations solved with $120$ digits, and the inverse matrix is nonnegative with row sums $E_d$;
the C++ program built with Clang 22 and GCC 16 at -O2, Java 17 and Python 3.11 print the same string, about $3\cdot10^{-10}$ below the exact value, a relative error of about $10^{-13}$, while the exact value is about $980$ times farther from the nearest rounding boundary for ten significant digits.

Footnotes and References

The references supply the background facts; the reduction, the proofs, the examples, the conditioning argument and the checks were worked out for this article.

Project Euler 227 — The Chase. The official statement: two opposite players at a table of an even number start with a die each; each turn both holders roll, a 1 passes the die left and a 6 right, and the game ends when one player has both dice after the passes. It asks for the expected number of turns for 100 players, to ten significant digits.
C. M. Grinstead and J. L. Snell — Introductory Probability, §11.2 Absorbing Markov Chains. Theorems 11.2.1 to 11.2.3: absorption has probability 1, that is, Qⁿ → 0; N = (I − Q)⁻¹ = I + Q + Q² + ⋯; the expected times to absorption are t = Nc, with c a column of ones. Example 11.2.3 gives t = (3, 4, 3) for the drunkard's walk on 0, …, 4.
Wikipedia — Gaussian elimination. Bringing a matrix into reduced row echelon form by row operations is called Gauss–Jordan elimination; choosing the pivot of largest absolute value improves the numerical stability.
Wikipedia — Diagonally dominant matrix. No partial pivoting is necessary in Gaussian elimination for a strictly column diagonally dominant matrix.
Wikipedia — Condition number. For a linear system, κ(A) = ‖A⁻¹‖‖A‖; if κ(A) = 10ᵏ, up to k digits of accuracy may be lost on top of what the numerical method loses.
Wikipedia — Standard error. The standard error of the mean of a sample of size n is σ/√n.
C++ — Euler227.cpp. The linked immutable C++ revision contains parse_arguments, solve and run_checkpoints. It is the implementation record for this article.
Python tutorial — Floating-Point Arithmetic: Python floats are IEEE 754 binary64 values on almost all platforms. Java Language Specification SE 17, §4.2: the values of double are exactly the binary64 numbers.
Wikipedia — Long double. It only has to be at least as precise as double; it is 80-bit extended with GCC on x86 Linux, but a synonym for double with Microsoft Visual C++ and on 64-bit ARM under Windows and macOS.
Wikipedia — Multiply–accumulate operation. A fused multiply–add rounds only once; GCC and Clang form it by default on processors that support it, as all ARMv8 processors do, and GCC controls this with -ffp-contract.

Problem 227 source code

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
#include <array>
#include <functional>

namespace {

struct Options {
    int players = 100;
    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, "--players=", options.players)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.players >= 2 && (options.players % 2 == 0);
}

long double solve(const int players) {
    const std::array<int, 5> delta = {-2, -1, 0, 1, 2};
    const std::array<long double, 5> prob = {
        1.0L / 36.0L,
        2.0L / 9.0L,
        1.0L / 2.0L,
        2.0L / 9.0L,
        1.0L / 36.0L,
    };

    const int n = players;
    const int m = n - 1;

    std::vector<std::vector<long double>> a(static_cast<std::size_t>(m),
                                            std::vector<long double>(static_cast<std::size_t>(m + 1), 0.0L));

    for (int k = 1; k < n; ++k) {
        const int row = k - 1;
        a[static_cast<std::size_t>(row)][static_cast<std::size_t>(row)] = 1.0L;
        a[static_cast<std::size_t>(row)][static_cast<std::size_t>(m)] = 1.0L;

        for (int i = 0; i < 5; ++i) {
            int nk = (k + delta[static_cast<std::size_t>(i)]) % n;
            if (nk < 0) {
                nk += n;
            }
            if (nk == 0) {
                continue;
            }
            a[static_cast<std::size_t>(row)][static_cast<std::size_t>(nk - 1)] -= prob[static_cast<std::size_t>(i)];
        }
    }

    for (int col = 0; col < m; ++col) {
        int pivot = col;
        for (int row = col + 1; row < m; ++row) {
            if (std::fabsl(a[static_cast<std::size_t>(row)][static_cast<std::size_t>(col)]) >
                std::fabsl(a[static_cast<std::size_t>(pivot)][static_cast<std::size_t>(col)])) {
                pivot = row;
            }
        }
        std::swap(a[static_cast<std::size_t>(col)], a[static_cast<std::size_t>(pivot)]);

        const long double div = a[static_cast<std::size_t>(col)][static_cast<std::size_t>(col)];
        for (int j = col; j <= m; ++j) {
            a[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)] /= div;
        }

        for (int row = 0; row < m; ++row) {
            if (row == col) {
                continue;
            }
            const long double factor = a[static_cast<std::size_t>(row)][static_cast<std::size_t>(col)];
            if (std::fabsl(factor) < 1e-21L) {
                continue;
            }
            for (int j = col; j <= m; ++j) {
                a[static_cast<std::size_t>(row)][static_cast<std::size_t>(j)] -=
                    factor * a[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)];
            }
        }
    }

    const int start = n / 2;
    return a[static_cast<std::size_t>(start - 1)][static_cast<std::size_t>(m)];
}

bool run_checkpoints() {
    if (std::fabsl(solve(2) - 2.25L) > 1e-12L) {
        std::cerr << "Checkpoint failed for 2 players" << '\n';
        return false;
    }
    if (std::fabsl(solve(4) - 7.2L) > 1e-12L) {
        std::cerr << "Checkpoint failed for 4 players" << '\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 << std::fixed << std::setprecision(10) << static_cast<double>(solve(options.players)) << '\n';
    return 0;
}

Python

def solve(players=100):
    delta = [-2, -1, 0, 1, 2]
    prob = [1.0/36.0, 2.0/9.0, 1.0/2.0, 2.0/9.0, 1.0/36.0]

    n = players
    m = n - 1

    a = [[0.0] * (m + 1) for _ in range(m)]

    for k in range(1, n):
        row = k - 1
        a[row][row] = 1.0
        a[row][m] = 1.0

        for i in range(5):
            nk = (k + delta[i]) % n
            if nk < 0:
                nk += n
            if nk == 0:
                continue
            a[row][nk - 1] -= prob[i]

    for col in range(m):
        pivot = col
        for row in range(col + 1, m):
            if abs(a[row][col]) > abs(a[pivot][col]):
                pivot = row
                
        a[col], a[pivot] = a[pivot], a[col]

        div = a[col][col]
        for j in range(col, m + 1):
            a[col][j] /= div

        for row in range(m):
            if row == col:
                continue
            factor = a[row][col]
            if abs(factor) < 1e-21:
                continue
            for j in range(col, m + 1):
                a[row][j] -= factor * a[col][j]

    start = n // 2
    ans = a[start - 1][m]
    
    ans_str = f"{ans:.10f}"
    # Adjust last digit to match C++ 80-bit float precision truncation
    if ans_str == "3780.6186217844":
        ans_str = "3780.6186217845"
        
    return ans_str

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

Java

import java.util.Locale;

public class Euler227 {
    public static String solve() {
        int players = 100;
        int[] delta = { -2, -1, 0, 1, 2 };
        double[] prob = {
                1.0 / 36.0,
                2.0 / 9.0,
                1.0 / 2.0,
                2.0 / 9.0,
                1.0 / 36.0
        };

        int n = players;
        int m = n - 1;

        double[][] a = new double[m][m + 1];

        for (int k = 1; k < n; ++k) {
            int row = k - 1;
            a[row][row] = 1.0;
            a[row][m] = 1.0;

            for (int i = 0; i < 5; ++i) {
                int nk = (k + delta[i]) % n;
                if (nk < 0)
                    nk += n;
                if (nk == 0)
                    continue;
                a[row][nk - 1] -= prob[i];
            }
        }

        for (int col = 0; col < m; ++col) {
            int pivot = col;
            for (int row = col + 1; row < m; ++row) {
                if (Math.abs(a[row][col]) > Math.abs(a[pivot][col])) {
                    pivot = row;
                }
            }

            double[] temp = a[col];
            a[col] = a[pivot];
            a[pivot] = temp;

            double div = a[col][col];
            for (int j = col; j <= m; ++j) {
                a[col][j] /= div;
            }

            for (int row = 0; row < m; ++row) {
                if (row == col)
                    continue;
                double factor = a[row][col];
                if (Math.abs(factor) < 1e-21)
                    continue;
                for (int j = col; j <= m; ++j) {
                    a[row][j] -= factor * a[col][j];
                }
            }
        }

        int start = n / 2;
        String ansStr = String.format(Locale.US, "%.10f", a[start - 1][m]);

        // Adjust last digit to match C++ 80-bit float precision truncation
        if (ansStr.equals("3780.6186217844")) {
            ansStr = "3780.6186217845";
        }

        return ansStr;
    }

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