Problem 740: Secret Santa
View on Project EulerProject Euler Problem 740 Solution
EulerSolve provides an optimized solution for Project Euler Problem 740, Secret Santa, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary There are \(n\) labeled participants and initially two slips carrying each label. The participants act in the order \(1,2,\dots,n-1\). When participant \(k\) acts, every remaining slip labeled \(k\) is forbidden to that participant, and then two admissible slips are drawn uniformly without replacement. The quantity \(q(n)\) is the probability that when the final participant is reached, at least one slip carrying label \(n\) is still present. A direct history-by-history simulation would branch far too much, because the exact identity of many remaining slips quickly stops mattering. The key idea is to replace the full multiset of labels by a symmetry-reduced state that remembers only the information relevant for future exclusions and for the survival of label \(n\). Mathematical Approach Let \(\mathcal{F}_k(\rho,\mu,\sigma_1,\sigma_2)\) be the success probability from stage \(k\), where the state variables mean: \(\rho\): how many slips labeled \(k\) are still present, so they are forbidden at the current stage; \(\mu\): how many slips labeled \(n\) are still present; \(\sigma_1\): how many future intermediate labels currently have exactly one self-labeled slip left; \(\sigma_2\): how many future intermediate labels currently have exactly two self-labeled slips left....
Detailed mathematical approach
Problem Summary
There are \(n\) labeled participants and initially two slips carrying each label. The participants act in the order \(1,2,\dots,n-1\). When participant \(k\) acts, every remaining slip labeled \(k\) is forbidden to that participant, and then two admissible slips are drawn uniformly without replacement. The quantity \(q(n)\) is the probability that when the final participant is reached, at least one slip carrying label \(n\) is still present.
A direct history-by-history simulation would branch far too much, because the exact identity of many remaining slips quickly stops mattering. The key idea is to replace the full multiset of labels by a symmetry-reduced state that remembers only the information relevant for future exclusions and for the survival of label \(n\).
Mathematical Approach
Let \(\mathcal{F}_k(\rho,\mu,\sigma_1,\sigma_2)\) be the success probability from stage \(k\), where the state variables mean:
\(\rho\): how many slips labeled \(k\) are still present, so they are forbidden at the current stage;
\(\mu\): how many slips labeled \(n\) are still present;
\(\sigma_1\): how many future intermediate labels currently have exactly one self-labeled slip left;
\(\sigma_2\): how many future intermediate labels currently have exactly two self-labeled slips left.
Step 1: Compress the Remaining Box by Symmetry
At stage \(k\), the total number of remaining slips is
$$T_k=2(n-k+1).$$
We do not need to remember every label separately. For any future intermediate participant, only the number \(0\), \(1\), or \(2\) of remaining self-labeled slips matters, because that number determines how many slips will be forbidden when that participant later becomes current. All labels in the same class are interchangeable.
The box therefore splits into four groups: the \(\rho\) forbidden slips of the current participant, the \(\mu\) slips of the final participant, the \(\sigma_1\) tracked slips belonging to class-1 future labels, and the \(2\sigma_2\) tracked slips belonging to class-2 future labels. Every other remaining slip can be merged into one generic pool. Its size is forced to be
$$\gamma=T_k-\rho-\mu-\sigma_1-2\sigma_2.$$
This generic pool contains slips whose exact labels no longer matter individually: slips of already processed participants, and also non-critical slips from future labels whose exclusion status is already encoded by \(\sigma_1\) and \(\sigma_2\).
Step 2: Describe One Stage as Two Draws Without Replacement
The current participant may draw from all slips except the \(\rho\) forbidden ones, so the admissible total is
$$A=T_k-\rho.$$
The four drawable categories have counts
$$c_0=\gamma,\qquad c_1=\mu,\qquad c_2=\sigma_1,\qquad c_3=2\sigma_2.$$
If the first draw uses category \(u\) and the second uses category \(v\), then the ordered probability is
$$\Pr(u,v)=\frac{c_u}{A}\cdot\frac{c_v^{(u)}}{A-1},$$
where \(c_v^{(u)}\) is the updated count after the first draw. Ordered pairs are necessary, because the second draw sees a different box.
Step 3: Update the Compressed Counts After a Draw
The effect of one draw depends only on its category:
$$\begin{aligned} \text{generic draw} &:\quad \gamma\to\gamma-1,\\ \text{final-label draw} &:\quad \mu\to\mu-1,\\ \text{class-1 draw} &:\quad \sigma_1\to\sigma_1-1,\\ \text{class-2 draw} &:\quad \sigma_2\to\sigma_2-1,\qquad \sigma_1\to\sigma_1+1. \end{aligned}$$
The last line reflects a future label that used to have two self-labeled slips and now has only one. After two draws we obtain updated values \(\mu',\sigma_1',\sigma_2'\), while \(\gamma\) can always be recomputed from the conservation formula.
Step 4: Average Over the Next Participant
After stage \(k\), the next current participant is \(k+1\). The compressed state does not record which future label is \(k+1\); it records only how many future labels lie in class \(0\), \(1\), or \(2\). By symmetry, \(k+1\) is uniformly distributed over those future labels.
Let
$$N_f=n-k-1,\qquad \tau_1=\sigma_1',\qquad \tau_2=\sigma_2',\qquad \tau_0=N_f-\tau_1-\tau_2.$$
If \(k=n-1\), there is no intermediate participant left, so the continuation is simply success or failure according to whether \(\mu'>0\).
Otherwise the continuation value is
$$\frac{\tau_0}{N_f}\mathcal{F}_{k+1}(0,\mu',\tau_1,\tau_2)+\frac{\tau_1}{N_f}\mathcal{F}_{k+1}(1,\mu',\tau_1-1,\tau_2)+\frac{\tau_2}{N_f}\mathcal{F}_{k+1}(2,\mu',\tau_1,\tau_2-1).$$
This is the central symmetry step: instead of tracking identities, we average over the three possible classes of the next participant.
Step 5: Boundary Condition and Initial State
When the process reaches stage \(n\), the only question left is whether at least one slip labeled \(n\) survived. Therefore
$$\mathcal{F}_n(\rho,\mu,\sigma_1,\sigma_2)= \begin{cases} 1,& \mu>0,\\ 0,& \mu=0. \end{cases}$$
At the start, participant \(1\) has both self-labeled slips still present, the final participant also has both of theirs, and every intermediate future participant starts in class \(2\). Hence
$$q(n)=\mathcal{F}_1(2,2,0,n-2).$$
Worked Example: \(n=3\)
The initial state is \(\mathcal{F}_1(2,2,0,1)\). Participant \(1\) cannot draw the two slips labeled \(1\), so the admissible box contains exactly four slips: two labeled \(2\) and two labeled \(3\).
If both draws remove label \(3\), the probability is
$$\frac{2}{4}\cdot\frac{1}{3}=\frac{1}{6},$$
and the process fails immediately because no slip labeled \(3\) remains.
If one draw removes label \(2\) and the other removes label \(3\), the total probability is
$$2\cdot\frac{2}{4}\cdot\frac{2}{3}=\frac{2}{3}.$$
Then participant \(2\) reaches stage \(2\) with one forbidden self-slip and one surviving slip labeled \(3\). Among the three admissible slips, only one carries label \(3\), so participant \(2\) leaves that slip untouched with probability \(1/3\).
If both draws remove label \(2\), the probability is again \(1/6\). Then participant \(2\) faces four admissible slips, two of which are labeled \(3\); the only failure is drawing both of them, so the continuation probability is \(1-1/6=5/6\).
Combining the three cases gives
$$q(3)=\frac{1}{6}\cdot 0+\frac{2}{3}\cdot\frac{1}{3}+\frac{1}{6}\cdot\frac{5}{6}=\frac{13}{36}=0.361111\ldots,$$
which matches the checkpoint used by the implementation.
How the Code Works
The C++, Python, and Java implementations all memoize the same five-parameter dynamic-programming state. The dense versions store values in a flattened table indexed by the stage and the four compressed counters, while the Python version uses a dictionary keyed by the same tuple. In every state, the implementation reconstructs the generic pool size from the conservation equation, computes the admissible total, and then loops over the \(4\times 4\) ordered category pairs for the two draws.
Pairs with zero count are skipped immediately. For each valid pair, the implementation applies the category updates twice, computes the continuation mixture over the three possible classes of the next participant, and adds the weighted contribution to the memoized answer. Because all symmetric labels are merged from the start, the program never needs to enumerate actual names or explicit draw histories.
Complexity Analysis
The stage index contributes \(O(n)\) possibilities, the two slip counters \(\rho\) and \(\mu\) each range over \(\{0,1,2\}\), and the two future-class counters range up to \(n\). Therefore the worst-case state space is \(O(n^4)\). Each state performs at most \(16\) ordered draw transitions and a constant-size continuation mixture, so the running time is \(O(n^4)\) with a moderate constant factor. The dense table versions use \(O(n^4)\) memory, while the memoized dictionary version stores only the reachable subset, though the same worst-case bound still applies.
Footnotes and References
- Problem page: https://projecteuler.net/problem=740
- Secret Santa: Wikipedia - Secret Santa
- Dynamic programming: Wikipedia - Dynamic programming
- Hypergeometric distribution: Wikipedia - Hypergeometric distribution
- Memoization: Wikipedia - Memoization
Problem 740 source code
C++
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
namespace {
struct Solver {
int n;
int dim;
std::vector<long double> memo;
std::vector<std::uint8_t> seen;
explicit Solver(const int n_) : n(n_), dim(n_ + 1) {
const std::size_t total = static_cast<std::size_t>(n + 1) * 3U * 3U *
static_cast<std::size_t>(dim) * static_cast<std::size_t>(dim);
memo.assign(total, 0.0L);
seen.assign(total, 0U);
}
std::size_t index(const int i, const int x, const int y, const int a1, const int a2) const {
std::size_t z = static_cast<std::size_t>(i);
z = z * 3U + static_cast<std::size_t>(x);
z = z * 3U + static_cast<std::size_t>(y);
z = z * static_cast<std::size_t>(dim) + static_cast<std::size_t>(a1);
z = z * static_cast<std::size_t>(dim) + static_cast<std::size_t>(a2);
return z;
}
static void apply_draw(const int cat, int& g, int& y, int& a1, int& a2) {
if (cat == 0) {
--g;
} else if (cat == 1) {
--y;
} else if (cat == 2) {
--a1;
} else {
--a2;
++a1;
}
}
long double solve(const int i, const int x, const int y, const int a1, const int a2) {
if (i == n) {
return (y > 0) ? 1.0L : 0.0L;
}
const std::size_t id = index(i, x, y, a1, a2);
if (seen[id] != 0U) {
return memo[id];
}
seen[id] = 1U;
const int future_pool = n - i - 1; // labels i+1..n-1
const int total_slips = 2 * (n - i + 1);
const int g0 = total_slips - x - y - a1 - 2 * a2;
const int avail0 = total_slips - x;
const auto count_for_cat = [](const int cat, const int g, const int yv, const int c1, const int c2) {
if (cat == 0) return g;
if (cat == 1) return yv;
if (cat == 2) return c1;
return 2 * c2;
};
long double ans = 0.0L;
for (int c1 = 0; c1 < 4; ++c1) {
const int cnt1 = count_for_cat(c1, g0, y, a1, a2);
if (cnt1 <= 0) {
continue;
}
const long double p1 = static_cast<long double>(cnt1) / static_cast<long double>(avail0);
int g1 = g0, y1 = y, a11 = a1, a21 = a2;
apply_draw(c1, g1, y1, a11, a21);
const int avail1 = avail0 - 1;
for (int c2 = 0; c2 < 4; ++c2) {
const int cnt2 = count_for_cat(c2, g1, y1, a11, a21);
if (cnt2 <= 0) {
continue;
}
const long double p2 = static_cast<long double>(cnt2) / static_cast<long double>(avail1);
int g2 = g1, y2 = y1, a12 = a11, a22 = a21;
apply_draw(c2, g2, y2, a12, a22);
long double cont = 0.0L;
if (i == n - 1) {
cont = (y2 > 0) ? 1.0L : 0.0L;
} else {
const int b1 = a12;
const int b2 = a22;
const int b0 = future_pool - b1 - b2;
if (b0 > 0) {
const long double ps = static_cast<long double>(b0) / static_cast<long double>(future_pool);
cont += ps * solve(i + 1, 0, y2, b1, b2);
}
if (b1 > 0) {
const long double ps = static_cast<long double>(b1) / static_cast<long double>(future_pool);
cont += ps * solve(i + 1, 1, y2, b1 - 1, b2);
}
if (b2 > 0) {
const long double ps = static_cast<long double>(b2) / static_cast<long double>(future_pool);
cont += ps * solve(i + 1, 2, y2, b1, b2 - 1);
}
}
ans += p1 * p2 * cont;
}
}
memo[id] = ans;
return ans;
}
};
long double q(const int n) {
Solver solver(n);
return solver.solve(1, 2, 2, 0, n - 2);
}
} // namespace
int main() {
assert(std::fabsl(q(3) - 0.3611111111L) < 5e-11L);
assert(std::fabsl(q(5) - 0.2476095994L) < 5e-11L);
std::cout << std::fixed << std::setprecision(10)
<< static_cast<double>(q(100)) << '\n';
return 0;
}
Python
def solve():
n = 100
dim = n + 1
memo = {}
def idx(i, x, y, a1, a2):
return (i, x, y, a1, a2)
def apply_draw(cat, g, y, a1, a2):
if cat == 0: return g-1, y, a1, a2
elif cat == 1: return g, y-1, a1, a2
elif cat == 2: return g, y, a1-1, a2
else: return g, y, a1+1, a2-1
def count_cat(cat, g, y, a1, a2):
if cat == 0: return g
elif cat == 1: return y
elif cat == 2: return a1
else: return 2*a2
def rec(i, x, y, a1, a2):
if i == n: return 1.0 if y > 0 else 0.0
key = (i, x, y, a1, a2)
if key in memo: return memo[key]
total_slips = 2*(n - i + 1)
g0 = total_slips - x - y - a1 - 2*a2
avail0 = total_slips - x
future = n - i - 1
ans = 0.0
for c1 in range(4):
cnt1 = count_cat(c1, g0, y, a1, a2)
if cnt1 <= 0: continue
p1 = cnt1 / avail0
g1, y1, a11, a21 = apply_draw(c1, g0, y, a1, a2)
avail1 = avail0 - 1
for c2 in range(4):
cnt2 = count_cat(c2, g1, y1, a11, a21)
if cnt2 <= 0: continue
p2 = cnt2 / avail1
g2, y2, a12, a22 = apply_draw(c2, g1, y1, a11, a21)
if i == n - 1:
cont = 1.0 if y2 > 0 else 0.0
else:
b1, b2 = a12, a22
b0 = future - b1 - b2
cont = 0.0
if b0 > 0: cont += (b0/future) * rec(i+1, 0, y2, b1, b2)
if b1 > 0: cont += (b1/future) * rec(i+1, 1, y2, b1-1, b2)
if b2 > 0: cont += (b2/future) * rec(i+1, 2, y2, b1, b2-1)
ans += p1 * p2 * cont
memo[key] = ans
return ans
return f'{rec(1, 2, 2, 0, n-2):.10f}'
if __name__ == '__main__':
print(solve())
Java
public class Euler740 {
static class Solver {
int n;
int dim;
double[] memo;
byte[] seen;
Solver(int n) {
this.n = n;
this.dim = n + 1;
int total = (n + 1) * 3 * 3 * dim * dim;
memo = new double[total];
seen = new byte[total];
}
int index(int i, int x, int y, int a1, int a2) {
return ((((i * 3 + x) * 3 + y) * dim + a1) * dim) + a2;
}
static int[] applyDraw(int cat, int g, int yv, int a1, int a2) {
if (cat == 0)
g--;
else if (cat == 1)
yv--;
else if (cat == 2)
a1--;
else {
a2--;
a1++;
}
return new int[] { g, yv, a1, a2 };
}
static int countForCat(int cat, int g, int yv, int c1, int c2) {
if (cat == 0)
return g;
if (cat == 1)
return yv;
if (cat == 2)
return c1;
return 2 * c2;
}
double solveDfs(int i, int x, int y, int a1, int a2) {
if (i == n) {
return (y > 0) ? 1.0 : 0.0;
}
int id = index(i, x, y, a1, a2);
if (seen[id] != 0) {
return memo[id];
}
seen[id] = 1;
int futurePool = n - i - 1;
int totalSlips = 2 * (n - i + 1);
int g0 = totalSlips - x - y - a1 - 2 * a2;
int avail0 = totalSlips - x;
double ans = 0.0;
for (int c1 = 0; c1 < 4; ++c1) {
int cnt1 = countForCat(c1, g0, y, a1, a2);
if (cnt1 <= 0)
continue;
double p1 = (double) cnt1 / avail0;
int[] state1 = applyDraw(c1, g0, y, a1, a2);
int g1 = state1[0], y1 = state1[1], a11 = state1[2], a21 = state1[3];
int avail1 = avail0 - 1;
for (int c2 = 0; c2 < 4; ++c2) {
int cnt2 = countForCat(c2, g1, y1, a11, a21);
if (cnt2 <= 0)
continue;
double p2 = (double) cnt2 / avail1;
int[] state2 = applyDraw(c2, g1, y1, a11, a21);
int g2 = state2[0], y2 = state2[1], a12 = state2[2], a22 = state2[3];
double cont = 0.0;
if (i == n - 1) {
cont = (y2 > 0) ? 1.0 : 0.0;
} else {
int b1 = a12;
int b2 = a22;
int b0 = futurePool - b1 - b2;
if (b0 > 0) {
double ps = (double) b0 / futurePool;
cont += ps * solveDfs(i + 1, 0, y2, b1, b2);
}
if (b1 > 0) {
double ps = (double) b1 / futurePool;
cont += ps * solveDfs(i + 1, 1, y2, b1 - 1, b2);
}
if (b2 > 0) {
double ps = (double) b2 / futurePool;
cont += ps * solveDfs(i + 1, 2, y2, b1, b2 - 1);
}
}
ans += p1 * p2 * cont;
}
}
memo[id] = ans;
return ans;
}
}
public static String solve() {
Solver solver = new Solver(100);
double ans = solver.solveDfs(1, 2, 2, 0, 98);
return String.format(java.util.Locale.US, "%.10f", ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}