Computes E(Delta | phi(k), n, m) exactly (no simulation) for n = 12345678, m = 12345, where f(k) = phi(k) (Euler's totient). S = sum_{k=1..n} f(k) S* = sum_{i=1..m} f(X_i) (X_i - X_{i-1}) for a random increasing m-tuple Delta = S - S* E(Delta) = sum_{k=1..n-m} f(k) * C(n-k, m) / C(n, m) with weights updated by recurrence w_{k+1} = w_k * (n-k-m)/(n-k). A rigorous tail bound truncates the sum so only the necessary initial part is computed, keeping the totient sieve small.
# Project Euler 756: Approximating a Sum
#
# Computes E(Delta | phi(k), n, m) exactly (no simulation) for
# n = 12345678, m = 12345, where f(k) = phi(k) (Euler's totient).
#
# S = sum_{k=1..n} f(k)
# S* = sum_{i=1..m} f(X_i) (X_i - X_{i-1}) for a random increasing m-tuple
# Delta = S - S*
#
# E(Delta) = sum_{k=1..n-m} f(k) * C(n-k, m) / C(n, m)
# with weights updated by recurrence w_{k+1} = w_k * (n-k-m)/(n-k).
#
# A rigorous tail bound truncates the sum so only the necessary initial
# part is computed, keeping the totient sieve small.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
# Euler's linear sieve: phi[k] = totient(k) for 0 <= k <= n.
function totients_up_to(n: i64) -> ptr<i32> {
let phi: ptr<i32> = calloc(n + 1, 4)
if phi == null {
return null
}
if n >= 1 {
phi[1] = 1
}
# primes list — generous cap (number of primes <= n is well under n).
let primes: ptr<i32> = calloc(n + 1, 4)
let mut count: i64 = 0
let mut i: i64 = 2
while i <= n {
if phi[i] == 0 {
# i is prime
primes[count] = (i as i32)
count = count + 1
phi[i] = ((i - 1) as i32)
}
let mut j: i64 = 0
while j < count {
let p: i64 = (primes[j] as i64)
let ip: i64 = i * p
if ip > n {
break
}
if i % p == 0 {
phi[ip] = (phi[i] as i64) * p
break
} else {
phi[ip] = (phi[i] as i64) * (p - 1)
}
j = j + 1
}
i = i + 1
}
free(primes)
return phi
}
# Find a safe truncation index K (<= n-m) such that the remaining tail
# contribution is provably < eps. Uses phi(k) <= k <= n and non-increasing
# weights: sum_{k>K} phi(k) w_k <= n * (n-m-K) * w_{K+1}.
function cutoff_index(n: i64, m: i64, eps: f64) -> i64 {
let limit: i64 = n - m
if limit <= 0 {
return 0
}
let mut w: f64 = ((n - m) as f64) / (n as f64)
let mut k: i64 = 1
while k <= limit {
let remaining: i64 = limit - k
let nk: i64 = n - k
let mut w_next: f64
if nk <= m {
w_next = 0.0
} else {
w_next = w * ((nk - m) as f64) / (nk as f64)
}
if (n as f64) * (remaining as f64) * w_next < eps {
return k
}
w = w_next
k = k + 1
}
return limit
}
function expected_error_for_phi(n: i64, m: i64) -> f64 {
let limit: i64 = n - m
if limit <= 0 {
return 0.0
}
let mut upto: i64 = cutoff_index(n, m, 5e-8)
if upto > limit {
upto = limit
}
let phi: ptr<i32> = totients_up_to(upto)
let mut w: f64 = ((n - m) as f64) / (n as f64)
let mut total: f64 = 0.0
let mut c: f64 = 0.0
let mut k: i64 = 1
while k <= upto {
let x: f64 = ((phi[k]) as f64) * w
# Kahan compensated summation
let y: f64 = x - c
let t: f64 = total + y
c = (t - total) - y
total = t
let nk: i64 = n - k
if nk <= m {
break
}
w = w * ((nk - m) as f64) / (nk as f64)
k = k + 1
}
free(phi)
return total
}
function main() -> i32 {
let n: i64 = 12345678
let m: i64 = 12345
printf("%.6f\n", expected_error_for_phi(n, m))
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; }
int32_t* totients_up_to_i64(int64_t n);
int64_t cutoff_index_i64_i64_f64(int64_t n, int64_t m, double eps);
double expected_error_for_phi_i64_i64(int64_t n, int64_t m);
int32_t main(void);
int32_t* totients_up_to_i64(int64_t n) {
int32_t* phi = (int32_t*)(calloc((n + 1), 4));
if (phi == NULL) {
return NULL;
}
if (n >= 1) {
phi[1] = 1;
}
int32_t* primes = (int32_t*)(calloc((n + 1), 4));
int64_t count = 0;
int64_t i = 2;
while (i <= n) {
if (phi[i] == 0) {
primes[count] = ((int32_t)(i));
count = (count + 1);
phi[i] = ((int32_t)((i - 1)));
}
int64_t j = 0;
while (j < count) {
int64_t p = ((int64_t)(primes[j]));
int64_t ip = (i * p);
if (ip > n) {
break;
}
if (FLOW_CHECKED_MOD((i), (p)) == 0) {
phi[ip] = (((int64_t)(phi[i])) * p);
break;
} else {
phi[ip] = (((int64_t)(phi[i])) * (p - 1));
}
j = (j + 1);
}
i = (i + 1);
}
free(primes);
return phi;
}
int64_t cutoff_index_i64_i64_f64(int64_t n, int64_t m, double eps) {
int64_t limit = (n - m);
if (limit <= 0) {
return 0;
}
double w = (((double)((n - m))) / ((double)(n)));
int64_t k = 1;
while (k <= limit) {
int64_t remaining = (limit - k);
int64_t nk = (n - k);
double w_next;
if (nk <= m) {
w_next = 0.0;
} else {
w_next = ((w * ((double)((nk - m)))) / ((double)(nk)));
}
if (((((double)(n)) * ((double)(remaining))) * w_next) < eps) {
return k;
}
w = w_next;
k = (k + 1);
}
return limit;
}
double expected_error_for_phi_i64_i64(int64_t n, int64_t m) {
int64_t limit = (n - m);
if (limit <= 0) {
return 0.0;
}
int64_t upto = cutoff_index_i64_i64_f64(n, m, 5e-8);
if (upto > limit) {
upto = limit;
}
int32_t* phi = (int32_t*)(totients_up_to_i64(upto));
double w = (((double)((n - m))) / ((double)(n)));
double total = 0.0;
double c = 0.0;
int64_t k = 1;
while (k <= upto) {
double x = (((double)(phi[k])) * w);
double y = (x - c);
double t = (total + y);
c = ((t - total) - y);
total = t;
int64_t nk = (n - k);
if (nk <= m) {
break;
}
w = ((w * ((double)((nk - m)))) / ((double)(nk)));
k = (k + 1);
}
free(phi);
return total;
}
int32_t main(void) {
int64_t n = 12345678;
int64_t m = 12345;
printf("%.6f\n", expected_error_for_phi_i64_i64(n, m));
return 0;
}