Problem 687: Shuffling Cards
View on Project EulerProject Euler Problem 687 Solution
EulerSolve provides an optimized solution for Project Euler Problem 687, Shuffling Cards, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary A uniformly shuffled standard deck has \(52\) cards, arranged as \(13\) ranks with \(4\) cards of each rank. A rank is called perfect if no two cards of that rank appear in adjacent positions in the shuffled deck. If \(Y\) denotes the number of perfect ranks, the problem asks for $$\Pr\!\left(Y\in\{2,3,5,7,11,13\}\right).$$ Brute force over all \(52!\) permutations is hopeless, so the implementations reorganize the question into two layers: first compute the probability that every rank in a fixed subset is perfect, then recover the full distribution of \(Y\) from those subset probabilities. Mathematical Approach The key observation is that the exact identities of the tracked ranks do not matter; only how many ranks are being tracked matters. That symmetry makes it possible to solve thirteen smaller probability problems and then combine them. Step 1: Fix a subset of ranks Choose a subset \(S\) of \(j\) ranks, where \(0\le j\le 13\), and define $$q_j=\Pr(\text{every rank in }S\text{ is perfect}).$$ Because every rank has the same multiplicity and every shuffle is uniform, \(q_j\) depends only on \(j\), not on the actual labels of the ranks in \(S\). If we can evaluate \(q_j\) for all \(j\), then we can rebuild the distribution of the random variable \(Y\)....
Detailed mathematical approach
Problem Summary
A uniformly shuffled standard deck has \(52\) cards, arranged as \(13\) ranks with \(4\) cards of each rank. A rank is called perfect if no two cards of that rank appear in adjacent positions in the shuffled deck. If \(Y\) denotes the number of perfect ranks, the problem asks for
$$\Pr\!\left(Y\in\{2,3,5,7,11,13\}\right).$$
Brute force over all \(52!\) permutations is hopeless, so the implementations reorganize the question into two layers: first compute the probability that every rank in a fixed subset is perfect, then recover the full distribution of \(Y\) from those subset probabilities.
Mathematical Approach
The key observation is that the exact identities of the tracked ranks do not matter; only how many ranks are being tracked matters. That symmetry makes it possible to solve thirteen smaller probability problems and then combine them.
Step 1: Fix a subset of ranks
Choose a subset \(S\) of \(j\) ranks, where \(0\le j\le 13\), and define
$$q_j=\Pr(\text{every rank in }S\text{ is perfect}).$$
Because every rank has the same multiplicity and every shuffle is uniform, \(q_j\) depends only on \(j\), not on the actual labels of the ranks in \(S\).
If we can evaluate \(q_j\) for all \(j\), then we can rebuild the distribution of the random variable \(Y\).
Step 2: Turn subset probabilities into binomial moments
For each fixed \(j\)-subset \(S\), let \(I_S\) be the indicator that all ranks in \(S\) are perfect. In any particular shuffle with exactly \(Y\) perfect ranks, the number of \(j\)-subsets that satisfy this is \(\binom{Y}{j}\). Therefore
$$\sum_{\lvert S\rvert=j} I_S=\binom{Y}{j}.$$
Taking expectations and using symmetry gives
$$B_j=\binom{13}{j}q_j=\mathbb{E}\!\left[\binom{Y}{j}\right].$$
Now write
$$D_k=\Pr(Y=k),\qquad 0\le k\le 13.$$
Then the \(j\)-th binomial moment satisfies
$$B_j=\sum_{k=j}^{13}\binom{k}{j}D_k.$$
This is an upper-triangular linear system, so once the \(B_j\) values are known, the exact probabilities \(D_k\) can be recovered by backward substitution.
Step 3: Reveal the deck sequentially and compress the state
To compute \(q_j\), imagine revealing the deck from left to right. Cards whose ranks are outside \(S\) are irrelevant except that they can separate two tracked cards and thereby prevent an adjacency violation.
At any point, the state can be summarized by:
$$u=\text{number of remaining cards from untracked ranks},$$
$$a_r=\text{number of tracked ranks with exactly }r\text{ unseen cards left},\qquad r=1,2,3,4,$$
and
$$s=\text{remaining multiplicity of the most recently drawn tracked rank},$$
with the convention \(s=0\) if the previous card was untracked or if the previous tracked rank has already been exhausted.
This compressed description is sufficient because ranks in the same class are exchangeable. The only rank that needs to be distinguished is the most recently drawn tracked rank, since drawing it again immediately would create an adjacent equal-rank pair and destroy perfection for that rank.
The initial state for a subset of size \(j\) is
$$\Psi_j(52-4j,0,0,0,j,0),$$
where \(\Psi_j\) denotes the probability of eventually finishing the reveal with every tracked rank still perfect.
Step 4: Write the transition probabilities
From a state \((u,a_1,a_2,a_3,a_4,s)\), the number of cards still unseen is
$$R=u+a_1+2a_2+3a_3+4a_4.$$
If \(R=0\), the tracked subset has survived the whole shuffle without any forbidden adjacency, so
$$\Psi_j(0,0,0,0,0,0)=1.$$
If an untracked card is drawn, this happens with probability \(u/R\), the count \(u\) drops by \(1\), and the active restriction disappears because an untracked card separates future tracked cards from the previously drawn tracked rank.
If a tracked rank currently belongs to class \(r\), then each such rank contributes \(r\) possible next cards. However, if \(s=r\), one of those ranks is the forbidden rank that cannot be repeated immediately, so the number of eligible ranks in class \(r\) is
$$a_r-\mathbf{1}_{\{s=r\}}.$$
Hence the probability of drawing a tracked card from class \(r\) is
$$\frac{r\left(a_r-\mathbf{1}_{\{s=r\}}\right)}{R}.$$
After such a draw, one tracked rank moves from class \(r\) to class \(r-1\). The new active value becomes \(r-1\), except that when \(r=1\) the rank is exhausted and the active value resets to \(0\).
Each state has at most five outgoing transitions: one for an untracked card and one for each multiplicity class \(r=1,2,3,4\). Memoization makes this recurrence practical because the same compressed state can be reached through many different reveal histories.
Step 5: Recover the exact distribution of \(Y\)
Once \(q_j=\Psi_j(52-4j,0,0,0,j,0)\) is known for every \(j\), we compute
$$B_j=\binom{13}{j}q_j.$$
The triangular relation
$$B_j=\sum_{k=j}^{13}\binom{k}{j}D_k$$
can then be inverted from \(k=13\) downward:
$$D_k=B_k-\sum_{t=k+1}^{13}\binom{t}{k}D_t.$$
Finally, the required probability is
$$\Pr\!\left(Y\in\{2,3,5,7,11,13\}\right)=D_2+D_3+D_5+D_7+D_{11}+D_{13}.$$
Worked Example: A single fixed rank
For \(j=1\), let us track one specific rank. It is perfect exactly when its four positions in the shuffled deck contain no adjacent pair.
If those positions are
$$1\le x_1<x_2<x_3<x_4\le 52,\qquad x_{i+1}\ge x_i+2,$$
then setting \(y_i=x_i-(i-1)\) transforms them into
$$1\le y_1<y_2<y_3<y_4\le 49.$$
So the number of non-adjacent \(4\)-position choices is \(\binom{49}{4}\), while the total number of \(4\)-position choices is \(\binom{52}{4}\). Therefore
$$q_1=\frac{\binom{49}{4}}{\binom{52}{4}}=\frac{4324}{5525}.$$
By linearity of expectation,
$$\mathbb{E}[Y]=13q_1=\frac{4324}{425}.$$
This matches the consistency check built into the implementations and confirms that the subset-probability viewpoint is aligned with the combinatorics of the shuffle.
How the Code Works
The C++, Python, and Java implementations all follow the same sequence. First they precompute the binomial coefficients \(\binom{n}{k}\) for \(0\le n,k\le 13\). Then, for each subset size \(j\), they evaluate the memoized probability recurrence starting from the initial state with \(j\) tracked ranks and \(52-4j\) untracked cards.
Those thirteen probabilities produce the values \(q_1,\dots,q_{13}\), and from them the implementations form the binomial moments \(B_j=\binom{13}{j}q_j\). The triangular system is then solved backwards to obtain the full distribution \(D_0,D_1,\dots,D_{13}\).
After the distribution is known, the implementation adds the six probabilities corresponding to prime values of \(Y\) and prints the result to ten decimal places. One implementation also checks the identity \(\mathbb{E}[Y]=4324/425\), which is exactly the worked example above written as a sanity test.
Complexity Analysis
For a fixed subset size \(j\), let \(S_j\) be the number of reachable compressed states \((u,a_1,a_2,a_3,a_4,s)\). Memoization ensures that each such state is evaluated once, and each evaluation considers at most five transitions. Therefore the running time for that \(j\) is \(O(S_j)\) and the memory usage is also \(O(S_j)\).
Across all subset sizes, the total cost is
$$O\!\left(\sum_{j=1}^{13} S_j\right)+O(13^2),$$
where the \(O(13^2)\) term is the final backward substitution. There is no simple closed form for \(S_j\), but the state space is tightly constrained by
$$u+a_1+2a_2+3a_3+4a_4\le 52$$
and by the fact that only \(13\) ranks exist. In practice the reachable-state count is small enough that the exact probability can be computed comfortably.
Footnotes and References
- Project Euler problem page: https://projecteuler.net/problem=687
- Binomial coefficient: Wikipedia - Binomial coefficient
- Linearity of expectation: Wikipedia - Expected value
- Dynamic programming: Wikipedia - Dynamic programming
- Binomial transform and inversion: Wikipedia - Binomial transform
Problem 687 source code
C++
#include <array>
#include <cassert>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <unordered_map>
#include <cmath>
#include <functional>
namespace {
using ld = long double;
struct ProbDP {
std::unordered_map<std::uint32_t, ld> memo;
static std::uint32_t encode(int f, int c1, int c2, int c3, int c4, int prev) {
const int p = (prev < 0) ? 0 : prev;
return static_cast<std::uint32_t>(f) |
(static_cast<std::uint32_t>(c1) << 6U) |
(static_cast<std::uint32_t>(c2) << 10U) |
(static_cast<std::uint32_t>(c3) << 14U) |
(static_cast<std::uint32_t>(c4) << 18U) |
(static_cast<std::uint32_t>(p) << 22U);
}
ld dfs(int f, int c1, int c2, int c3, int c4, int prev) {
const int rem = f + c1 + 2 * c2 + 3 * c3 + 4 * c4;
if (rem == 0) {
return 1.0L;
}
const std::uint32_t key = encode(f, c1, c2, c3, c4, prev);
const auto it = memo.find(key);
if (it != memo.end()) {
return it->second;
}
ld out = 0.0L;
if (f > 0) {
out += (static_cast<ld>(f) / static_cast<ld>(rem)) *
dfs(f - 1, c1, c2, c3, c4, -1);
}
auto add_rank = [&](int r, int count, int nc1, int nc2, int nc3, int nc4) {
if (count == 0) {
return;
}
int avail = count;
if (prev == r) {
--avail;
}
if (avail <= 0) {
return;
}
const int np = (r - 1 == 0) ? -1 : (r - 1);
const ld prob = static_cast<ld>(avail * r) / static_cast<ld>(rem);
out += prob * dfs(f, nc1, nc2, nc3, nc4, np);
};
add_rank(1, c1, c1 - 1, c2, c3, c4);
add_rank(2, c2, c1 + 1, c2 - 1, c3, c4);
add_rank(3, c3, c1, c2 + 1, c3 - 1, c4);
add_rank(4, c4, c1, c2, c3 + 1, c4 - 1);
memo.emplace(key, out);
return out;
}
ld prob_all_perfect_for_subset_size(int j) {
memo.clear();
memo.reserve(1 << 18);
return dfs(52 - 4 * j, 0, 0, 0, j, -1);
}
};
} // namespace
int main() {
std::array<std::array<int, 14>, 14> C{};
for (int n = 0; n <= 13; ++n) {
C[static_cast<std::size_t>(n)][0] = 1;
C[static_cast<std::size_t>(n)][static_cast<std::size_t>(n)] = 1;
for (int k = 1; k < n; ++k) {
C[static_cast<std::size_t>(n)][static_cast<std::size_t>(k)] =
C[static_cast<std::size_t>(n - 1)][static_cast<std::size_t>(k - 1)] +
C[static_cast<std::size_t>(n - 1)][static_cast<std::size_t>(k)];
}
}
ProbDP solver;
std::array<ld, 14> p{};
p[0] = 1.0L;
for (int j = 1; j <= 13; ++j) {
p[static_cast<std::size_t>(j)] = solver.prob_all_perfect_for_subset_size(j);
}
std::array<ld, 14> M{};
for (int j = 0; j <= 13; ++j) {
M[static_cast<std::size_t>(j)] =
static_cast<ld>(C[13][j]) * p[static_cast<std::size_t>(j)];
}
std::array<ld, 14> P{};
for (int k = 13; k >= 0; --k) {
ld v = M[static_cast<std::size_t>(k)];
for (int t = k + 1; t <= 13; ++t) {
v -= static_cast<ld>(C[static_cast<std::size_t>(t)][static_cast<std::size_t>(k)]) *
P[static_cast<std::size_t>(t)];
}
P[static_cast<std::size_t>(k)] = v;
}
ld expect = 0.0L;
for (int k = 0; k <= 13; ++k) {
expect += static_cast<ld>(k) * P[static_cast<std::size_t>(k)];
}
assert(std::fabsl(expect - (4324.0L / 425.0L)) < 1e-12L);
const std::array<int, 6> primes{{2, 3, 5, 7, 11, 13}};
ld ans = 0.0L;
for (int pval : primes) {
ans += P[static_cast<std::size_t>(pval)];
}
std::cout << std::fixed << std::setprecision(10) << static_cast<double>(ans) << "\n";
return 0;
}
Python
def solve():
memo = {}
def dfs(f, c1, c2, c3, c4, prev):
rem = f + c1 + 2 * c2 + 3 * c3 + 4 * c4
if rem == 0:
return 1.0
key = (f, c1, c2, c3, c4, prev)
if key in memo:
return memo[key]
out = 0.0
if f > 0:
out += (f / rem) * dfs(f - 1, c1, c2, c3, c4, -1)
def add_rank(r, count, nc1, nc2, nc3, nc4):
nonlocal out
if count == 0:
return
avail = count
if prev == r:
avail -= 1
if avail <= 0:
return
np_val = -1 if r == 1 else r - 1
prob = (avail * r) / rem
out += prob * dfs(f, nc1, nc2, nc3, nc4, np_val)
add_rank(1, c1, c1 - 1, c2, c3, c4)
add_rank(2, c2, c1 + 1, c2 - 1, c3, c4)
add_rank(3, c3, c1, c2 + 1, c3 - 1, c4)
add_rank(4, c4, c1, c2, c3 + 1, c4 - 1)
memo[key] = out
return out
def prob_all_perfect(j):
memo.clear()
return dfs(52 - 4 * j, 0, 0, 0, j, -1)
C = [[0] * 14 for _ in range(14)]
for n in range(14):
C[n][0] = 1
C[n][n] = 1
for k in range(1, n):
C[n][k] = C[n - 1][k - 1] + C[n - 1][k]
p = [0.0] * 14
p[0] = 1.0
for j in range(1, 14):
p[j] = prob_all_perfect(j)
M = [0.0] * 14
for j in range(14):
M[j] = C[13][j] * p[j]
P = [0.0] * 14
for k in range(13, -1, -1):
v = M[k]
for t in range(k + 1, 14):
v -= C[t][k] * P[t]
P[k] = v
primes = [2, 3, 5, 7, 11, 13]
ans = 0.0
for pval in primes:
ans += P[pval]
return f"{ans:.10f}"
if __name__ == '__main__':
print(solve())
Java
import java.util.HashMap;
import java.util.Map;
import java.util.Locale;
public class Euler687 {
static class ProbDP {
Map<Integer, Double> memo = new HashMap<>();
static int encode(int f, int c1, int c2, int c3, int c4, int prev) {
int p = (prev < 0) ? 0 : prev;
return f | (c1 << 6) | (c2 << 10) | (c3 << 14) | (c4 << 18) | (p << 22);
}
double dfs(int f, int c1, int c2, int c3, int c4, int prev) {
int rem = f + c1 + 2 * c2 + 3 * c3 + 4 * c4;
if (rem == 0)
return 1.0;
int key = encode(f, c1, c2, c3, c4, prev);
if (memo.containsKey(key)) {
return memo.get(key);
}
double out = 0.0;
if (f > 0) {
out += ((double) f / rem) * dfs(f - 1, c1, c2, c3, c4, -1);
}
if (c1 > 0) {
int avail = c1;
if (prev == 1)
avail--;
if (avail > 0) {
double prob = (double) (avail * 1) / rem;
out += prob * dfs(f, c1 - 1, c2, c3, c4, -1);
}
}
if (c2 > 0) {
int avail = c2;
if (prev == 2)
avail--;
if (avail > 0) {
double prob = (double) (avail * 2) / rem;
out += prob * dfs(f, c1 + 1, c2 - 1, c3, c4, 1);
}
}
if (c3 > 0) {
int avail = c3;
if (prev == 3)
avail--;
if (avail > 0) {
double prob = (double) (avail * 3) / rem;
out += prob * dfs(f, c1, c2 + 1, c3 - 1, c4, 2);
}
}
if (c4 > 0) {
int avail = c4;
if (prev == 4)
avail--;
if (avail > 0) {
double prob = (double) (avail * 4) / rem;
out += prob * dfs(f, c1, c2, c3 + 1, c4 - 1, 3);
}
}
memo.put(key, out);
return out;
}
double probAllPerfect(int j) {
memo.clear();
return dfs(52 - 4 * j, 0, 0, 0, j, -1);
}
}
public static String solve() {
int[][] C = new int[14][14];
for (int n = 0; n <= 13; ++n) {
C[n][0] = 1;
C[n][n] = 1;
for (int k = 1; k < n; ++k) {
C[n][k] = C[n - 1][k - 1] + C[n - 1][k];
}
}
ProbDP solver = new ProbDP();
double[] p = new double[14];
p[0] = 1.0;
for (int j = 1; j <= 13; ++j) {
p[j] = solver.probAllPerfect(j);
}
double[] M = new double[14];
for (int j = 0; j <= 13; ++j) {
M[j] = C[13][j] * p[j];
}
double[] P = new double[14];
for (int k = 13; k >= 0; --k) {
double v = M[k];
for (int t = k + 1; t <= 13; ++t) {
v -= C[t][k] * P[t];
}
P[k] = v;
}
int[] primes = { 2, 3, 5, 7, 11, 13 };
double ans = 0.0;
for (int pval : primes) {
ans += P[pval];
}
return String.format(Locale.US, "%.10f", ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}