Expected stopping time for drawing cards without replacement from 10 decks (540 cards) until every suit, rank, and deck design has appeared. Uses inclusion-exclusion over missing suits (a), ranks (b), and deck designs (c). For each combination the allowed card count is M(a,b,c) = (4-a)*(13-b)*(10-c) + 2*(10-c) and the coefficient is (-1)^(a+b+c+1) * C(4,a)*C(13,b)*C(10,c). E[T] = sum_M coeff[M] * S(M, N) where S(M,N) = sum_{k=0}^{N-1} C(M,k)/C(N,k), computed via the recurrence r_0 = 1, r_k = r_{k-1} * (M-k+1)/(N-k+1). Uses i128 fixed-point arithmetic (scale 10^24) for precision.
# Project Euler 796: A Grand Shuffle
#
# Expected stopping time for drawing cards without replacement from 10 decks
# (540 cards) until every suit, rank, and deck design has appeared.
#
# Uses inclusion-exclusion over missing suits (a), ranks (b), and deck
# designs (c). For each combination the allowed card count is
# M(a,b,c) = (4-a)*(13-b)*(10-c) + 2*(10-c)
# and the coefficient is (-1)^(a+b+c+1) * C(4,a)*C(13,b)*C(10,c).
#
# E[T] = sum_M coeff[M] * S(M, N)
# where S(M,N) = sum_{k=0}^{N-1} C(M,k)/C(N,k), computed via the recurrence
# r_0 = 1, r_k = r_{k-1} * (M-k+1)/(N-k+1).
#
# Uses i128 fixed-point arithmetic (scale 10^24) for precision.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
function pow10(n: i32) -> i128 {
let mut result: i128 = 1
for i in 0..n {
result = result * (10 as i128)
}
return result
}
function binom(n: i32, k: i32) -> i64 {
if k < 0 || k > n {
return 0
}
let mut kk: i32 = k
if kk > n - kk {
kk = n - kk
}
let mut result: i64 = 1
for i in 0..kk {
result = result * ((n - i) as i64) / ((i + 1) as i64)
}
return result
}
function main() -> i32 {
let N: i64 = 540
let max_M: i32 = 540
let SCALE: i128 = pow10(24)
# Compute coefficients via inclusion-exclusion.
let coeffs: ptr<i64> = calloc((max_M + 1) as i64, 8)
if coeffs == null {
return 1
}
for i in 0..((max_M + 1) as i64) {
coeffs[i] = 0
}
for a in 0..5 {
for b in 0..14 {
for c in 0..11 {
if a == 0 && b == 0 && c == 0 {
continue
}
let mut coef: i64 = binom(4, a) * binom(13, b) * binom(10, c)
if (a + b + c) % 2 == 0 {
coef = -coef
}
let decks_left: i32 = 10 - c
let M: i32 = (4 - a) * (13 - b) * decks_left + 2 * decks_left
coeffs[M] = coeffs[M] + coef
}
}
}
# Compute S(M, N) for each M with nonzero coefficient, in fixed-point.
let mut total_scaled: i128 = 0
for M in 0..(max_M + 1) {
if coeffs[M] == 0 {
continue
}
let mut r: i128 = SCALE
let mut s_val: i128 = SCALE
let mut max_k: i64 = N - 1
if (M as i64) < max_k {
max_k = M as i64
}
for k in 1..(max_k + 1) {
r = (r * ((M as i64 - k + 1) as i128)) / ((N - k + 1) as i128)
s_val = s_val + r
}
total_scaled = total_scaled + (coeffs[M] as i128) * s_val
}
# Round to 8 decimal places and output.
# rounded = (total_scaled + 5*10^15) / 10^16
let half_round: i128 = pow10(15) * (5 as i128)
let round_div: i128 = pow10(16)
let ten_p_8: i128 = pow10(8)
let rounded: i128 = (total_scaled + half_round) / round_div
let int_part: i64 = (rounded / ten_p_8) as i64
let frac_part: i64 = (rounded % ten_p_8) as i64
printf("%lld.%08lld\n", int_part, frac_part)
free(coeffs)
return 0
}
Generated C
#include <stdint.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/* Flow runtime helpers */
typedef struct flow_temp_node { struct flow_temp_node* next; } flow_temp_node;
static flow_temp_node* flow_temp_head = NULL;
static int flow_temp_atexit_set = 0;
__attribute__((unused)) static void flow_temp_free_all(void) {
while (flow_temp_head) {
flow_temp_node* n = flow_temp_head;
flow_temp_head = n->next;
free(n);
}
}
__attribute__((unused)) static void* flow_temp_alloc(size_t nbytes) {
flow_temp_node* node = (flow_temp_node*)malloc(sizeof(flow_temp_node) + nbytes);
if (!node) return NULL;
node->next = flow_temp_head;
flow_temp_head = node;
if (!flow_temp_atexit_set) {
flow_temp_atexit_set = 1;
atexit(flow_temp_free_all);
}
return (void*)(node + 1);
}
#ifndef FLOW_DIAG
#define FLOW_DIAG(msg) fprintf(stderr, "%s", (msg))
#endif
#ifndef FLOW_LOG
#define FLOW_LOG(fmt, ...) printf(fmt, __VA_ARGS__)
#endif
#ifndef FLOW_LOG_EMPTY
#define FLOW_LOG_EMPTY(fmt) printf(fmt)
#endif
static char* flow_strcat(const char* a, const char* b) {
size_t la = strlen(a ? a : ""), lb = strlen(b ? b : "");
char* r = (char*)flow_temp_alloc(la + lb + 1);
if (!r) return NULL;
if (la) memcpy(r, a, la);
if (lb) memcpy(r + la, b, lb);
r[la + lb] = '\0';
return r;
}
#define __flow_in_arr(arr, val) __extension__ ({ \
int _found = 0; \
size_t _n = sizeof(arr)/sizeof((arr)[0]); \
for (size_t _i = 0; _i < _n; _i++) { \
if ((arr)[_i] == (val)) { _found = 1; break; } \
} _found; })
/* Unified fault handler (MISRA #279) — override with -DFLOW_FAULT_HANDLER=fn */
#ifndef FLOW_FAULT_HANDLER
__attribute__((unused)) static inline void flow_fault_handler(const char* msg) {
fprintf(stderr, "flow: %s\n", msg ? msg : "fault");
abort();
#if defined(__GNUC__) || defined(__clang__)
__builtin_unreachable();
#endif
}
#else
#define flow_fault_handler FLOW_FAULT_HANDLER
#endif
#define flow_div_by_zero_handler() flow_fault_handler("division by zero")
#define flow_shift_ub_handler() flow_fault_handler("invalid shift (amount out of range or left-shift of negative)")
#ifndef FLOW_CHECKED_DIV
#define FLOW_CHECKED_DIV(L, R) (((R) != 0) ? ((L) / (R)) : (flow_div_by_zero_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_MOD
#define FLOW_CHECKED_MOD(L, R) (((R) != 0) ? ((L) % (R)) : (flow_div_by_zero_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_SHL
#define FLOW_CHECKED_SHL(L, R) ((((R) >= 0) && ((unsigned long long)(R) < (sizeof(L) * 8ull)) && ((L) >= 0)) ? ((L) << (R)) : (flow_shift_ub_handler(), (L) * 0))
#endif
#ifndef FLOW_CHECKED_SHR
#define FLOW_CHECKED_SHR(L, R) ((((R) >= 0) && ((unsigned long long)(R) < (sizeof(L) * 8ull))) ? ((L) >> (R)) : (flow_shift_ub_handler(), (L) * 0))
#endif
#include <math.h>
void* _ui_state = NULL;
static inline float i32_to_f32(int32_t v) { return (float)v; }
/* Host stub for @gpu kernels (device codegen replaces this). */
static inline int32_t gpu_thread_id(void) { return 0; }
__int128 pow10_i32(int32_t n);
int64_t binom_i32_i32(int32_t n, int32_t k);
int32_t main(void);
__int128 pow10_i32(int32_t n) {
__int128 result = 1;
int32_t __flow_step_1 = 1;
for (int32_t i = 0; (0 <= n) ? i < n : i > n; i += (0 <= n) ? 1 : -1) {
result = (result * ((__int128)(10)));
}
return result;
}
int64_t binom_i32_i32(int32_t n, int32_t k) {
if ((k < 0 || k > n)) {
return 0;
}
int32_t kk = k;
if (kk > (n - kk)) {
kk = (n - kk);
}
int64_t result = 1;
int32_t __flow_step_2 = 1;
for (int32_t i = 0; (0 <= kk) ? i < kk : i > kk; i += (0 <= kk) ? 1 : -1) {
result = FLOW_CHECKED_DIV(((result * ((int64_t)((n - i))))), (((int64_t)((i + 1)))));
}
return result;
}
int32_t main(void) {
int64_t N = 540;
int32_t max_M = 540;
__int128 SCALE = pow10_i32(24);
int64_t* coeffs = (int64_t*)(calloc(((int64_t)((max_M + 1))), 8));
if (coeffs == NULL) {
return 1;
}
int32_t __flow_step_3 = 1;
for (int32_t i = 0; (0 <= ((int64_t)((max_M + 1)))) ? i < ((int64_t)((max_M + 1))) : i > ((int64_t)((max_M + 1))); i += (0 <= ((int64_t)((max_M + 1)))) ? 1 : -1) {
coeffs[i] = 0;
}
int32_t __flow_step_4 = 1;
for (int32_t a = 0; (0 <= 5) ? a < 5 : a > 5; a += (0 <= 5) ? 1 : -1) {
int32_t __flow_step_5 = 1;
for (int32_t b = 0; (0 <= 14) ? b < 14 : b > 14; b += (0 <= 14) ? 1 : -1) {
int32_t __flow_step_6 = 1;
for (int32_t c = 0; (0 <= 11) ? c < 11 : c > 11; c += (0 <= 11) ? 1 : -1) {
if (((a == 0 && b == 0) && c == 0)) {
continue;
}
int64_t coef = ((binom_i32_i32(4, a) * binom_i32_i32(13, b)) * binom_i32_i32(10, c));
if (FLOW_CHECKED_MOD((((a + b) + c)), (2)) == 0) {
coef = (-coef);
}
int32_t decks_left = (10 - c);
int32_t M = ((((4 - a) * (13 - b)) * decks_left) + (2 * decks_left));
coeffs[M] = (coeffs[M] + coef);
}
}
}
__int128 total_scaled = 0;
int32_t __flow_step_7 = 1;
for (int32_t M = 0; (0 <= (max_M + 1)) ? M < (max_M + 1) : M > (max_M + 1); M += (0 <= (max_M + 1)) ? 1 : -1) {
if (coeffs[M] == 0) {
continue;
}
__int128 r = SCALE;
__int128 s_val = SCALE;
int64_t max_k = (N - 1);
if (((int64_t)(M)) < max_k) {
max_k = ((int64_t)(M));
}
int32_t __flow_step_8 = 1;
for (int32_t k = 1; (1 <= (max_k + 1)) ? k < (max_k + 1) : k > (max_k + 1); k += (1 <= (max_k + 1)) ? 1 : -1) {
r = FLOW_CHECKED_DIV(((r * ((__int128)(((((int64_t)(M)) - k) + 1))))), (((__int128)(((N - k) + 1)))));
s_val = (s_val + r);
}
total_scaled = (total_scaled + (((__int128)(coeffs[M])) * s_val));
}
__int128 half_round = (pow10_i32(15) * ((__int128)(5)));
__int128 round_div = pow10_i32(16);
__int128 ten_p_8 = pow10_i32(8);
__int128 rounded = FLOW_CHECKED_DIV(((total_scaled + half_round)), (round_div));
int64_t int_part = ((int64_t)(FLOW_CHECKED_DIV((rounded), (ten_p_8))));
int64_t frac_part = ((int64_t)(FLOW_CHECKED_MOD((rounded), (ten_p_8))));
printf("%lld.%08lld\n", int_part, frac_part);
free(coeffs);
return 0;
}