Problem 499: St. Petersburg Lottery
View on Project EulerProject Euler Problem 499 Solution
EulerSolve provides an optimized solution for Project Euler Problem 499, St. Petersburg Lottery, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Each round costs \(t\) units. The payout is \(2^n\) with probability \(2^{-(n+1)}\) for \(n \ge 0\). Starting from bankroll \(B\), we want the probability of being able to keep playing forever, which is the same as never letting the bankroll fall below the amount needed for the next ticket. The bankroll can in principle grow without bound, so the state space is infinite. The implementation therefore replaces the raw infinite recursion by a finite residue model together with a self-similar matrix equation. Mathematical Approach The numerical method follows the same probability law in all three implementations. Its key idea is that after shifting the bankroll by \(t-1\), the process becomes translation-invariant in blocks of size \(t-1\). Step 1: Shift the State and Truncate the Tiny Tail Let $$d=t-1,\qquad c=B-d.$$ Here \(c\) is the bankroll measured above the minimum playable threshold. One round changes \(c\) by $$c' = c + 2^n - t.$$ If \(c' \le 0\), the player can no longer buy another ticket, so that path is ruined. The exact distribution has infinitely many outcomes \(n=0,1,2,\dots\), but the implementation keeps only $$0 \le n \lt H,\qquad H=50.$$ The omitted probability mass is $$\sum_{n=H}^{\infty} 2^{-(n+1)} = 2^{-50},$$ which is far below the requested \(7\)-decimal accuracy....
Detailed mathematical approach
Problem Summary
Each round costs \(t\) units. The payout is \(2^n\) with probability \(2^{-(n+1)}\) for \(n \ge 0\). Starting from bankroll \(B\), we want the probability of being able to keep playing forever, which is the same as never letting the bankroll fall below the amount needed for the next ticket. The bankroll can in principle grow without bound, so the state space is infinite. The implementation therefore replaces the raw infinite recursion by a finite residue model together with a self-similar matrix equation.
Mathematical Approach
The numerical method follows the same probability law in all three implementations. Its key idea is that after shifting the bankroll by \(t-1\), the process becomes translation-invariant in blocks of size \(t-1\).
Step 1: Shift the State and Truncate the Tiny Tail
Let
$$d=t-1,\qquad c=B-d.$$
Here \(c\) is the bankroll measured above the minimum playable threshold. One round changes \(c\) by
$$c' = c + 2^n - t.$$
If \(c' \le 0\), the player can no longer buy another ticket, so that path is ruined.
The exact distribution has infinitely many outcomes \(n=0,1,2,\dots\), but the implementation keeps only
$$0 \le n \lt H,\qquad H=50.$$
The omitted probability mass is
$$\sum_{n=H}^{\infty} 2^{-(n+1)} = 2^{-50},$$
which is far below the requested \(7\)-decimal accuracy.
Step 2: Decompose the State into Layer and Residue
Every playable shifted bankroll \(c \ge 1\) can be written uniquely as
$$c = k d + r + 1,\qquad k \ge 0,\qquad 0 \le r \lt d.$$
The integer \(k\) is the layer, and \(r\) is the residue inside that layer. This decomposition is useful because
$$k d + r + 1 + 2^n - t = (k-1)d + (r + 2^n).$$
So increasing the starting layer by one simply shifts all future layers by one. The detailed geometry inside each layer stays the same.
Define
$$u_k(r)=\Pr(\text{eventual ruin} \mid c = k d + r + 1).$$
For fixed \(t\), each \(u_k\) is a vector of length \(d\).
Step 3: Build the Transfer Matrix for Higher Layers
Because the process repeats the same pattern from one layer to the next, there is a single \(d \times d\) matrix \(X\) such that
$$u_{k+1} = X u_k,\qquad u_k = X^k u_0.$$
To derive \(X\), start from layer \(1\), where \(c = d + r + 1\). After outcome \(n\),
$$c' = d + r + 1 + 2^n - t = r + 2^n.$$
Write this new state as
$$c' = j_n(r)\,d + \rho_n(r) + 1,$$
with
$$j_n(r)=\left\lfloor\frac{r+2^n-1}{d}\right\rfloor,\qquad \rho_n(r)=(r+2^n-1)\bmod d.$$
If \(j_n(r)=0\), the process lands in the boundary layer. If \(j_n(r)\ge 1\), it lands \(j_n(r)\) layers above the boundary. Therefore the \(r\)-th row of \(X\) satisfies
$$X_{r,*}=\sum_{j_n(r)=0} 2^{-(n+1)} e_{\rho_n(r)}^{\mathsf T}+\sum_{j_n(r)\ge 1} 2^{-(n+1)}\bigl(X^{j_n(r)}\bigr)_{\rho_n(r),*},$$
where \(e_{\rho}^{\mathsf T}\) is the row vector with a \(1\) in column \(\rho\) and zeros elsewhere.
This is a matrix fixed-point equation. The implementations iterate it numerically until successive matrices differ by less than \(10^{-18}\), or until the safety iteration cap is reached.
Step 4: Solve the Boundary Layer Separately
Now consider the lowest playable layer, where
$$c = r + 1,\qquad 0 \le r \lt d.$$
After outcome \(n\),
$$c' = r + 1 + 2^n - t = r + 2^n - d.$$
If \(c' \le 0\), ruin happens immediately. Otherwise write
$$c' = j'_n(r)\,d + \rho'_n(r) + 1,$$
where
$$j'_n(r)=\left\lfloor\frac{c'-1}{d}\right\rfloor,\qquad \rho'_n(r)=(c'-1)\bmod d.$$
Conditioning on the first round gives
$$u_0(r)=q_r+\sum_{j'_n(r)=0}2^{-(n+1)}u_0\!\left(\rho'_n(r)\right)+\sum_{j'_n(r)\ge 1}2^{-(n+1)}\bigl[X^{j'_n(r)}u_0\bigr]_{\rho'_n(r)},$$
where \(q_r\) is the total probability of immediate ruin from residue \(r\).
All boundary-layer unknowns can therefore be grouped into a finite linear system
$$(I-C)u_0=q,$$
which the implementations solve by Gaussian elimination with pivoting.
Step 5: Recover the Required Probability
For the given starting bankroll \(B\), first convert it to shifted capital
$$c_0=B-d.$$
If \(c_0 \le 0\), the answer is \(0\) because the first ticket cannot even be bought. Otherwise write
$$c_0 = k d + r + 1.$$
The ruin probability is then
$$\bigl[X^k u_0\bigr]_r,$$
so the desired survival probability is
$$\boxed{P_t(B)=1-\bigl[X^k u_0\bigr]_r.}$$
Step 6: Worked Example for \(t=2\)
When \(t=2\), we have \(d=1\), so there is only one residue class. The transfer matrix becomes a single scalar \(x\), and the boundary ruin vector becomes a single scalar \(u\).
From layer \(1\), outcome \(n=0\) sends the process to the boundary layer, while \(n\ge 1\) jumps to layer \(2^n-1\). Therefore
$$x=\frac12+\sum_{n=1}^{H-1}2^{-(n+1)}x^{2^n-1}.$$
For the boundary layer, \(n=0\) is immediate ruin and \(n\ge 1\) sends the process to layer \(2^n-2\), so
$$u=\frac12+\sum_{n=1}^{H-1}2^{-(n+1)}x^{2^n-2}u.$$
Thus
$$P_2(2)=1-u,\qquad P_2(5)=1-x^3u.$$
The numerical solution gives
$$P_2(2)\approx 0.2522,\qquad P_2(5)\approx 0.6873,$$
which matches the validation checkpoints used by the implementations.
How the Code Works
The C++, Python, and Java implementations all follow the same numerical plan. First they precompute the retained payouts \(2^n\) and their probabilities \(2^{-(n+1)}\) for \(0 \le n \lt 50\). Next they enumerate every transition from one reference layer and record only the resulting residue together with the layer jump that must be applied later.
To evaluate many powers of the same transfer matrix efficiently, the implementation builds \(X\), \(X^2\), \(X^4\), and so on, and reconstructs each required \(X^j\) from the binary expansion of \(j\). After the fixed-point iteration for \(X\) stabilizes, the boundary layer is folded into one finite linear system and solved by Gaussian elimination. Finally the initial bankroll is decomposed into layer and residue, the needed matrix power is applied to the boundary vector, and the result is subtracted from \(1\). The Python entry point is only a thin wrapper around the same compiled numerical core.
Complexity Analysis
Let \(d=t-1\). Suppose the reference-layer and boundary transitions involve \(E\) distinct positive layer jumps, and let \(J_{\max}\) be the largest of them. One fixed-point iteration needs the squaring chain up to
$$L=\left\lfloor\log_2 J_{\max}\right\rfloor+1,$$
so its dense-matrix cost is \(O((L+E)d^3)\). Solving the boundary linear system costs \(O(d^3)\). Once the transfer matrix is known, evaluating the answer for a starting bankroll in layer \(k\) needs \(O(d^3 \log k)\) time by binary exponentiation. Memory usage is \(O((L+E)d^2)\). For the actual target \(t=15\), this is entirely manageable because \(d=14\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=499
- St. Petersburg paradox: Wikipedia — St. Petersburg paradox
- Gambler's ruin: Wikipedia — Gambler's ruin
- Numerical note: truncating the payout law at \(H=50\) leaves tail probability \(2^{-50}\), which is negligible at \(7\)-decimal precision.
- Implementation note: the fixed-point iteration targets entrywise change below \(10^{-18}\), and the boundary problem is solved as a finite linear system.
Problem 499 source code
C++
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <unordered_map>
#include <vector>
using namespace std;
struct Matrix {
int n = 0;
vector<long double> a;
Matrix() = default;
explicit Matrix(int n_) : n(n_), a(static_cast<size_t>(n_) * n_, 0.0L) {}
long double& operator()(int r, int c) {
return a[static_cast<size_t>(r) * n + c];
}
long double operator()(int r, int c) const {
return a[static_cast<size_t>(r) * n + c];
}
};
static Matrix identity_matrix(int n) {
Matrix I(n);
for (int i = 0; i < n; ++i) {
I(i, i) = 1.0L;
}
return I;
}
static Matrix multiply(const Matrix& A, const Matrix& B) {
int n = A.n;
Matrix C(n);
for (int i = 0; i < n; ++i) {
for (int k = 0; k < n; ++k) {
long double aik = A(i, k);
if (aik == 0.0L) {
continue;
}
const long double* brow = &B.a[static_cast<size_t>(k) * n];
long double* crow = &C.a[static_cast<size_t>(i) * n];
for (int j = 0; j < n; ++j) {
crow[j] += aik * brow[j];
}
}
}
return C;
}
static vector<long double> multiply_vec(const Matrix& A, const vector<long double>& v) {
int n = A.n;
vector<long double> res(n, 0.0L);
for (int i = 0; i < n; ++i) {
long double sum = 0.0L;
const long double* row = &A.a[static_cast<size_t>(i) * n];
for (int j = 0; j < n; ++j) {
sum += row[j] * v[j];
}
res[i] = sum;
}
return res;
}
static vector<Matrix> build_pow2(const Matrix& R, int bits) {
vector<Matrix> pow2;
pow2.reserve(bits);
pow2.push_back(R);
for (int i = 1; i < bits; ++i) {
pow2.push_back(multiply(pow2.back(), pow2.back()));
}
return pow2;
}
static Matrix pow_from_bits(const vector<Matrix>& pow2, const vector<int>& bits, int n) {
Matrix res = identity_matrix(n);
for (int bit : bits) {
res = multiply(res, pow2[bit]);
}
return res;
}
static bool solve_linear_system(Matrix A, vector<long double> b, vector<long double>& x) {
int n = A.n;
x.assign(n, 0.0L);
for (int i = 0; i < n; ++i) {
int pivot = i;
long double max_abs = fabsl(A(i, i));
for (int r = i + 1; r < n; ++r) {
long double val = fabsl(A(r, i));
if (val > max_abs) {
max_abs = val;
pivot = r;
}
}
if (max_abs < 1e-30L) {
return false;
}
if (pivot != i) {
for (int c = i; c < n; ++c) {
swap(A(i, c), A(pivot, c));
}
swap(b[i], b[pivot]);
}
long double diag = A(i, i);
for (int c = i; c < n; ++c) {
A(i, c) /= diag;
}
b[i] /= diag;
for (int r = 0; r < n; ++r) {
if (r == i) {
continue;
}
long double factor = A(r, i);
if (factor == 0.0L) {
continue;
}
for (int c = i; c < n; ++c) {
A(r, c) -= factor * A(i, c);
}
b[r] -= factor * b[i];
}
}
x = b;
return true;
}
struct Transition {
int from = 0;
int to = 0;
int exp_index = 0;
long double prob = 0.0L;
};
struct ExpInfo {
long long exp = 0;
vector<int> bits;
};
class LotterySolver {
public:
LotterySolver(long long m_val, int max_n)
: m(m_val), d(static_cast<int>(m_val - 1)), N(max_n) {
build_transitions();
}
long double probability_no_ruin(long long s) {
if (s < m) {
return 0.0L;
}
Matrix R = compute_R();
vector<long double> y0 = compute_y0(R);
long long c = s - d;
long long k = (c - 1) / d;
int r = static_cast<int>((c - 1) % d);
if (k == 0) {
return 1.0L - y0[r];
}
Matrix Rk = matrix_power(R, k);
vector<long double> yk = multiply_vec(Rk, y0);
return 1.0L - yk[r];
}
private:
long long m;
int d;
int N;
vector<long long> pow2n;
vector<long double> probs;
Matrix A_minus;
vector<Transition> transitions;
vector<ExpInfo> exp_info;
void build_transitions() {
pow2n.resize(N);
probs.resize(N);
for (int n = 0; n < N; ++n) {
pow2n[n] = 1LL << n;
probs[n] = ldexp(1.0L, -(n + 1));
}
A_minus = Matrix(d);
vector<long long> exps;
transitions.clear();
for (int n = 0; n < N; ++n) {
long long delta = pow2n[n] - m;
long double prob = probs[n];
for (int r = 0; r < d; ++r) {
long long c = d + r + 1;
long long cprime = c + delta;
long long kprime = (cprime - 1) / d;
long long i = kprime - 1;
int rprime = static_cast<int>((cprime - 1) % d);
if (i == -1) {
A_minus(r, rprime) += prob;
} else {
long long exp = i + 1;
exps.push_back(exp);
transitions.push_back({r, rprime, 0, prob});
}
}
}
sort(exps.begin(), exps.end());
exps.erase(unique(exps.begin(), exps.end()), exps.end());
unordered_map<long long, int> exp_index;
exp_info.clear();
exp_info.reserve(exps.size());
for (size_t i = 0; i < exps.size(); ++i) {
exp_index[exps[i]] = static_cast<int>(i);
ExpInfo info;
info.exp = exps[i];
long long e = exps[i];
int bit = 0;
while (e > 0) {
if (e & 1LL) {
info.bits.push_back(bit);
}
e >>= 1LL;
++bit;
}
exp_info.push_back(info);
}
size_t idx = 0;
for (int n = 0; n < N; ++n) {
long long delta = pow2n[n] - m;
long double prob = probs[n];
for (int r = 0; r < d; ++r) {
long long c = d + r + 1;
long long cprime = c + delta;
long long kprime = (cprime - 1) / d;
long long i = kprime - 1;
if (i == -1) {
continue;
}
long long exp = i + 1;
transitions[idx].exp_index = exp_index[exp];
transitions[idx].prob = prob;
++idx;
}
}
}
Matrix matrix_power(const Matrix& R, long long exp) const {
if (exp == 0) {
return identity_matrix(d);
}
int max_bit = 0;
long long temp = exp;
while (temp > 0) {
++max_bit;
temp >>= 1LL;
}
vector<Matrix> pow2 = build_pow2(R, max_bit);
vector<int> bits;
long long e = exp;
int bit = 0;
while (e > 0) {
if (e & 1LL) {
bits.push_back(bit);
}
e >>= 1LL;
++bit;
}
return pow_from_bits(pow2, bits, d);
}
Matrix compute_R() {
Matrix R = A_minus;
const long double tol = 1e-18L;
const int max_iter = 3000;
int max_bit = 0;
for (const auto& info : exp_info) {
if (!info.bits.empty()) {
max_bit = max(max_bit, info.bits.back() + 1);
}
}
if (max_bit == 0) {
max_bit = 1;
}
for (int iter = 0; iter < max_iter; ++iter) {
vector<Matrix> pow2 = build_pow2(R, max_bit);
vector<Matrix> Rexp(exp_info.size(), Matrix(d));
for (size_t i = 0; i < exp_info.size(); ++i) {
Rexp[i] = pow_from_bits(pow2, exp_info[i].bits, d);
}
Matrix Rnew = A_minus;
for (const auto& tr : transitions) {
const Matrix& mat = Rexp[tr.exp_index];
const long double* src_row = &mat.a[static_cast<size_t>(tr.to) * d];
long double* dst_row = &Rnew.a[static_cast<size_t>(tr.from) * d];
for (int j = 0; j < d; ++j) {
dst_row[j] += tr.prob * src_row[j];
}
}
long double diff = 0.0L;
for (size_t i = 0; i < R.a.size(); ++i) {
diff = max(diff, fabsl(Rnew.a[i] - R.a[i]));
}
R = Rnew;
if (diff < tol) {
break;
}
if (iter == max_iter - 1) {
cerr << "Warning: R iteration did not fully converge (diff=" << setprecision(6) << diff << ")\n";
}
}
return R;
}
vector<long double> compute_y0(const Matrix& R) {
vector<long long> boundary_exps;
boundary_exps.reserve(static_cast<size_t>(2 * N));
for (int n = 0; n < N; ++n) {
long long delta = pow2n[n] - m;
for (int r = 0; r < d; ++r) {
long long c = r + 1;
long long cprime = c + delta;
if (cprime <= 0) {
continue;
}
long long kprime = (cprime - 1) / d;
if (kprime >= 1) {
boundary_exps.push_back(kprime);
}
}
}
sort(boundary_exps.begin(), boundary_exps.end());
boundary_exps.erase(unique(boundary_exps.begin(), boundary_exps.end()), boundary_exps.end());
unordered_map<long long, int> exp_index;
vector<vector<int>> exp_bits(boundary_exps.size());
int max_bit = 0;
for (size_t i = 0; i < boundary_exps.size(); ++i) {
exp_index[boundary_exps[i]] = static_cast<int>(i);
long long e = boundary_exps[i];
int bit = 0;
while (e > 0) {
if (e & 1LL) {
exp_bits[i].push_back(bit);
}
e >>= 1LL;
++bit;
}
if (!exp_bits[i].empty()) {
max_bit = max(max_bit, exp_bits[i].back() + 1);
}
}
if (max_bit == 0) {
max_bit = 1;
}
vector<Matrix> pow2 = build_pow2(R, max_bit);
vector<Matrix> Rexp(boundary_exps.size(), Matrix(d));
for (size_t i = 0; i < boundary_exps.size(); ++i) {
Rexp[i] = pow_from_bits(pow2, exp_bits[i], d);
}
Matrix M(d);
vector<long double> b(d, 0.0L);
for (int n = 0; n < N; ++n) {
long long delta = pow2n[n] - m;
long double prob = probs[n];
for (int r = 0; r < d; ++r) {
long long c = r + 1;
long long cprime = c + delta;
if (cprime <= 0) {
b[r] += prob;
continue;
}
long long kprime = (cprime - 1) / d;
int rprime = static_cast<int>((cprime - 1) % d);
if (kprime == 0) {
M(r, rprime) += prob;
} else {
int idx = exp_index[kprime];
const Matrix& mat = Rexp[idx];
const long double* src_row = &mat.a[static_cast<size_t>(rprime) * d];
long double* dst_row = &M.a[static_cast<size_t>(r) * d];
for (int j = 0; j < d; ++j) {
dst_row[j] += prob * src_row[j];
}
}
}
}
Matrix A = identity_matrix(d);
for (int i = 0; i < d; ++i) {
for (int j = 0; j < d; ++j) {
A(i, j) -= M(i, j);
}
}
vector<long double> y0;
if (!solve_linear_system(A, b, y0)) {
cerr << "Failed to solve boundary system.\n";
y0.assign(d, 1.0L);
}
return y0;
}
};
static void run_check(long long m, long long s, long double expected, long double tol) {
const int N = 50;
LotterySolver solver(m, N);
long double res = solver.probability_no_ruin(s);
long double diff = fabsl(res - expected);
cout << "p_" << m << "(" << s << ") = " << fixed << setprecision(7) << res;
if (diff <= tol) {
cout << " [PASS]";
} else {
cout << " [FAIL]";
}
cout << "\n";
}
int main() {
cout << "--- Validation Checkpoints ---\n";
run_check(2, 2, 0.2522L, 5e-4L);
run_check(2, 5, 0.6873L, 5e-4L);
run_check(6, 10000, 0.9952L, 5e-4L);
cout << "\n--- Final Solution ---\n";
const long long m = 15;
const long long s = 1000000000LL;
const int N = 50;
LotterySolver solver(m, N);
long double ans = solver.probability_no_ruin(s);
cout << fixed << setprecision(7) << ans << "\n";
cout << "Answer: " << fixed << setprecision(7) << ans << "\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
import java.util.*;
public class Euler499 {
static double[][] matMul(double[][] A, double[][] B) {
int n = A.length;
double[][] C = new double[n][n];
for (int i = 0; i < n; i++)
for (int k = 0; k < n; k++) {
if (A[i][k] == 0)
continue;
for (int j = 0; j < n; j++)
C[i][j] += A[i][k] * B[k][j];
}
return C;
}
static double[] matVec(double[][] A, double[] v) {
int n = A.length;
double[] r = new double[n];
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
r[i] += A[i][j] * v[j];
return r;
}
static double[][] matPow(double[][] A, long exp) {
int n = A.length;
double[][] res = new double[n][n];
for (int i = 0; i < n; i++)
res[i][i] = 1;
while (exp > 0) {
if ((exp & 1) == 1)
res = matMul(res, A);
A = matMul(A, A);
exp >>= 1;
}
return res;
}
static boolean solveLinear(double[][] A, double[] b, double[] x) {
int n = A.length;
for (int i = 0; i < n; i++) {
int piv = i;
for (int r = i + 1; r < n; r++)
if (Math.abs(A[r][i]) > Math.abs(A[piv][i]))
piv = r;
if (Math.abs(A[piv][i]) < 1e-30)
return false;
double[] tmp = A[i];
A[i] = A[piv];
A[piv] = tmp;
double t = b[i];
b[i] = b[piv];
b[piv] = t;
double d = A[i][i];
for (int c = i; c < n; c++)
A[i][c] /= d;
b[i] /= d;
for (int r = 0; r < n; r++) {
if (r == i)
continue;
double f = A[r][i];
for (int c = i; c < n; c++)
A[r][c] -= f * A[i][c];
b[r] -= f * b[i];
}
}
System.arraycopy(b, 0, x, 0, n);
return true;
}
static double probabilityNoRuin(long m, long s) {
if (s < m)
return 0;
int d = (int) (m - 1), N = 50;
long[] pow2n = new long[N];
double[] probs = new double[N];
for (int n = 0; n < N; n++) {
pow2n[n] = 1L << n;
probs[n] = Math.pow(2, -(n + 1));
}
// Build A_minus and transitions for R iteration
double[][] Aminus = new double[d][d];
List<int[]> trIdx = new ArrayList<>();
List<Double> trProb = new ArrayList<>();
Map<Long, Integer> expIndex = new TreeMap<>();
List<Long> exps = new ArrayList<>();
for (int n = 0; n < N; n++) {
long delta = pow2n[n] - m;
for (int r = 0; r < d; r++) {
long c = d + r + 1, cp = c + delta, kp = (cp - 1) / d;
int rp = (int) ((cp - 1) % d);
if (kp - 1 == -1)
Aminus[r][rp] += probs[n];
else {
long e = kp;
if (!expIndex.containsKey(e)) {
expIndex.put(e, exps.size());
exps.add(e);
}
trIdx.add(new int[] { r, rp, expIndex.get(e) });
trProb.add(probs[n]);
}
}
}
// Iterate R to fixed point
double[][] R = new double[d][d];
for (int i = 0; i < d; i++)
System.arraycopy(Aminus[i], 0, R[i], 0, d);
for (int iter = 0; iter < 500; iter++) {
double[][][] Rpows = new double[exps.size()][][];
for (int i = 0; i < exps.size(); i++)
Rpows[i] = matPow(R, exps.get(i));
double[][] Rn = new double[d][d];
for (int i = 0; i < d; i++)
System.arraycopy(Aminus[i], 0, Rn[i], 0, d);
for (int t = 0; t < trIdx.size(); t++) {
int[] ti = trIdx.get(t);
double p = trProb.get(t);
double[][] mp = Rpows[ti[2]];
for (int j = 0; j < d; j++)
Rn[ti[0]][j] += p * mp[ti[1]][j];
}
double diff = 0;
for (int i = 0; i < d; i++)
for (int j = 0; j < d; j++)
diff = Math.max(diff, Math.abs(Rn[i][j] - R[i][j]));
R = Rn;
if (diff < 1e-18)
break;
}
// Compute y0 boundary
List<Long> bexps = new ArrayList<>();
for (int n = 0; n < N; n++) {
long delta = pow2n[n] - m;
for (int r = 0; r < d; r++) {
long cp = r + 1 + delta;
if (cp <= 0)
continue;
long kp = (cp - 1) / d;
if (kp >= 1)
bexps.add(kp);
}
}
bexps = new ArrayList<>(new TreeSet<>(bexps));
double[][][] bRpows = new double[bexps.size()][][];
Map<Long, Integer> bExpIdx = new HashMap<>();
for (int i = 0; i < bexps.size(); i++) {
bExpIdx.put(bexps.get(i), i);
bRpows[i] = matPow(R, bexps.get(i));
}
double[][] M = new double[d][d];
double[] bb = new double[d];
for (int n = 0; n < N; n++) {
long delta = pow2n[n] - m;
for (int r = 0; r < d; r++) {
long cp = r + 1 + delta;
if (cp <= 0) {
bb[r] += probs[n];
continue;
}
long kp = (cp - 1) / d;
int rp = (int) ((cp - 1) % d);
if (kp == 0)
M[r][rp] += probs[n];
else {
int idx = bExpIdx.get(kp);
for (int j = 0; j < d; j++)
M[r][j] += probs[n] * bRpows[idx][rp][j];
}
}
}
double[][] A = new double[d][d];
for (int i = 0; i < d; i++) {
for (int j = 0; j < d; j++)
A[i][j] = -M[i][j];
A[i][i] += 1;
}
double[] y0 = new double[d];
solveLinear(A, bb, y0);
long c2 = s - d;
long k = (c2 - 1) / d;
int r2 = (int) ((c2 - 1) % d);
if (k == 0)
return 1.0 - y0[r2];
double[][] Rk = matPow(R, k);
double[] yk = matVec(Rk, y0);
return 1.0 - yk[r2];
}
public static void main(String[] args) {
double ans = probabilityNoRuin(15, 1000000000L);
System.out.printf("%.7f%n", ans);
}
}