# Project Euler 595
# Incremental Random Sort
# S(52) expected shuffles, rounded to 8 decimals (f64 DP).
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
function floor(x: f64) -> f64
}
function expected_shuffles(n: i64) -> f64 {
let fact: ptr<f64> = calloc(n + 1, 8)
if fact == null { return -1.0 }
fact[0] = 1.0
fact[1] = 1.0
let mut i: i64 = 2
while i <= n {
fact[i] = fact[i - 1] * (i as f64)
i = i + 1
}
# a[m*(n) + r] succession counts; store rows packed with stride n
let a: ptr<f64> = calloc((n + 1) * n, 8)
if a == null { free(fact); return -1.0 }
a[1 * n + 0] = 1.0
let mut m: i64 = 2
while m <= n {
let mut r: i64 = 0
while r < m {
# C(m-1, r)
let mut c1: f64 = 1.0
let mut t: i64 = 0
let mut kk: i64 = r
if kk > m - 1 - kk { kk = m - 1 - kk }
t = 1
while t <= kk {
c1 = c1 * ((m - 1 - kk + t) as f64) / (t as f64)
t = t + 1
}
let mut s: f64 = 0.0
let mut j: i64 = 0
while j < m - r {
let k: i64 = m - r - j
# C(m-1-r, j) * k!
let mut cj: f64 = 1.0
let mut jj: i64 = j
if jj > m - 1 - r - jj { jj = m - 1 - r - jj }
t = 1
while t <= jj {
cj = cj * ((m - 1 - r - jj + t) as f64) / (t as f64)
t = t + 1
}
let term: f64 = cj * fact[k]
if (j & 1) != 0 {
s = s - term
} else {
s = s + term
}
j = j + 1
}
a[m * n + r] = c1 * s
r = r + 1
}
m = m + 1
}
let T: ptr<f64> = calloc(n + 1, 8)
if T == null { free(a); free(fact); return -1.0 }
T[1] = 0.0
m = 2
while m <= n {
let denom: f64 = fact[m] - a[m * n + 0]
let mut num: f64 = fact[m]
let mut r: i64 = 1
while r < m {
num = num + a[m * n + r] * T[m - r]
r = r + 1
}
T[m] = num / denom
m = m + 1
}
let mut Sn: f64 = 0.0
let mut r: i64 = 0
while r < n {
Sn = Sn + (a[n * n + r] / fact[n]) * T[n - r]
r = r + 1
}
free(T)
free(a)
free(fact)
return Sn
}
function main() -> i32 {
let s: f64 = expected_shuffles(52)
# round half-up to 8 decimals
let scaled: f64 = s * 100000000.0
let mut q: i64 = floor(scaled) as i64
let frac: f64 = scaled - (q as f64)
if frac >= 0.5 {
q = q + 1
}
let ip: i64 = q / 100000000
let fp: i64 = q % 100000000
printf("%lld.%08lld\n", ip, fp)
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; }
double expected_shuffles_i64(int64_t n);
int32_t main(void);
double expected_shuffles_i64(int64_t n) {
double* fact = (double*)(calloc((n + 1), 8));
if (fact == NULL) {
return (-1.0);
}
fact[0] = 1.0;
fact[1] = 1.0;
int64_t i = 2;
while (i <= n) {
fact[i] = (fact[(i - 1)] * ((double)(i)));
i = (i + 1);
}
double* a = (double*)(calloc(((n + 1) * n), 8));
if (a == NULL) {
free(fact);
return (-1.0);
}
a[((1 * n) + 0)] = 1.0;
int64_t m = 2;
while (m <= n) {
int64_t r = 0;
while (r < m) {
double c1 = 1.0;
int64_t t = 0;
int64_t kk = r;
if (kk > ((m - 1) - kk)) {
kk = ((m - 1) - kk);
}
t = 1;
while (t <= kk) {
c1 = ((c1 * ((double)((((m - 1) - kk) + t)))) / ((double)(t)));
t = (t + 1);
}
double s = 0.0;
int64_t j = 0;
while (j < (m - r)) {
int64_t k = ((m - r) - j);
double cj = 1.0;
int64_t jj = j;
if (jj > (((m - 1) - r) - jj)) {
jj = (((m - 1) - r) - jj);
}
t = 1;
while (t <= jj) {
cj = ((cj * ((double)(((((m - 1) - r) - jj) + t)))) / ((double)(t)));
t = (t + 1);
}
double term = (cj * fact[k]);
if ((j & 1) != 0) {
s = (s - term);
} else {
s = (s + term);
}
j = (j + 1);
}
a[((m * n) + r)] = (c1 * s);
r = (r + 1);
}
m = (m + 1);
}
double* T = (double*)(calloc((n + 1), 8));
if (T == NULL) {
free(a);
free(fact);
return (-1.0);
}
T[1] = 0.0;
m = 2;
while (m <= n) {
double denom = (fact[m] - a[((m * n) + 0)]);
double num = fact[m];
int64_t r = 1;
while (r < m) {
num = (num + (a[((m * n) + r)] * T[(m - r)]));
r = (r + 1);
}
T[m] = (num / denom);
m = (m + 1);
}
double Sn = 0.0;
int64_t r = 0;
while (r < n) {
Sn = (Sn + ((a[((n * n) + r)] / fact[n]) * T[(n - r)]));
r = (r + 1);
}
free(T);
free(a);
free(fact);
return Sn;
}
int32_t main(void) {
double s = expected_shuffles_i64(52);
double scaled = (s * 100000000.0);
int64_t q = ((int64_t)(floor(scaled)));
double frac = (scaled - ((double)(q)));
if (frac >= 0.5) {
q = (q + 1);
}
int64_t ip = FLOW_CHECKED_DIV((q), (100000000));
int64_t fp = FLOW_CHECKED_MOD((q), (100000000));
printf("%lld.%08lld\n", ip, fp);
return 0;
}