# Project Euler 786
# Billiard -- B(10^9).
# Uses floor_sum, Mertens function with Du Jiao sieve, and Mobius inversion.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
function malloc(n: i64) -> ptr<void>
}
const HT_SIZE: i64 = 2097152
const HT_MASK: i64 = 2097151
let mut g_pref: ptr<i64> = null
let mut g_limit: i64 = 0
let mut g_ht_m_key: ptr<i64> = null
let mut g_ht_m_val: ptr<i64> = null
let mut g_ht_f_key: ptr<i64> = null
let mut g_ht_f_val: ptr<i64> = null
function floor_sum(n0: i64, m0: i64, a0: i64, b0: i64) -> i64 {
let mut res: i64 = 0
let mut n: i64 = n0
let mut m: i64 = m0
let mut a: i64 = a0
let mut b: i64 = b0
while true {
if a >= m {
res = res + (n - 1) * n * (a / m) / 2
a = a % m
}
if b >= m {
res = res + n * (b / m)
b = b % m
}
let y_max: i64 = a * n + b
if y_max < m { break }
b = y_max % m
n = y_max / m
let tmp: i64 = m
m = a
a = tmp
}
return res
}
function icbrt(n: i64) -> i64 {
if n <= 0 { return 0 }
let mut bits: i64 = 0
let mut t: i64 = n
while t != 0 {
bits = bits + 1
t = t >> 1
}
let mut x: i64 = 1 << ((bits + 2) / 3)
while true {
let y: i64 = (2 * x + n / (x * x)) / 3
if y >= x { break }
x = y
}
while (x + 1) * (x + 1) * (x + 1) <= n {
x = x + 1
}
while x * x * x > n {
x = x - 1
}
return x
}
function ht_get(ht_key: ptr<i64>, ht_val: ptr<i64>, key: i64, out: ptr<i64>) -> i64 {
let mut h: i64 = (key * 2654435761) & HT_MASK
while ht_key[h] != 0 {
if ht_key[h] == key {
out[0] = ht_val[h]
return 1
}
h = (h + 1) & HT_MASK
}
return 0
}
function ht_put(ht_key: ptr<i64>, ht_val: ptr<i64>, key: i64, val: i64) -> void {
let mut h: i64 = (key * 2654435761) & HT_MASK
while ht_key[h] != 0 {
if ht_key[h] == key {
ht_val[h] = val
return
}
h = (h + 1) & HT_MASK
}
ht_key[h] = key
ht_val[h] = val
}
function M_func(n: i64) -> i64 {
if n <= 0 { return 0 }
if n <= g_limit { return g_pref[n] }
let out: ptr<i64> = malloc(8) as ptr<i64>
if ht_get(g_ht_m_key, g_ht_m_val, n, out) == 1 {
let v: i64 = out[0]
free(out)
return v
}
free(out)
let mut res: i64 = 1
let mut l: i64 = 2
while l <= n {
let q: i64 = n / l
let r: i64 = n / q
res = res - (r - l + 1) * M_func(q)
l = r + 1
}
ht_put(g_ht_m_key, g_ht_m_val, n, res)
return res
}
function F_func(n: i64) -> i64 {
if n <= 0 { return 0 }
let out: ptr<i64> = malloc(8) as ptr<i64>
if ht_get(g_ht_f_key, g_ht_f_val, n, out) == 1 {
let v: i64 = out[0]
free(out)
return v
}
free(out)
let res: i64 = M_func(n) + F_func(n / 3)
ht_put(g_ht_f_key, g_ht_f_val, n, res)
return res
}
function count_points_nonprimitive(M: i64) -> i64 {
if M < 28 { return 0 }
let n: i64 = (M - 18) / 10
let b: i64 = M - 10 * n
let total: i64 = floor_sum(n, 18, 10, b)
let n3: i64 = n / 3
let b3: i64 = M - 30 * n3
let total3: i64 = floor_sum(n3, 18, 30, b3)
return total - total3
}
function count_points_primitive(M: i64) -> i64 {
let max_d: i64 = M / 28
if max_d <= 0 { return 0 }
let limit: i64 = icbrt(max_d * max_d) + 64
# Linear sieve for mu
let mu: ptr<i64> = calloc(limit + 1, 8) as ptr<i64>
let is_comp: ptr<i8> = calloc(limit + 1, 1) as ptr<i8>
let primes: ptr<i64> = malloc((limit + 1) * 8) as ptr<i64>
let mut pc: i64 = 0
mu[1] = 1
for i in 2..(limit + 1) {
if is_comp[i] == 0 {
primes[pc] = i
pc = pc + 1
mu[i] = -1
}
let mut j: i64 = 0
while j < pc {
let p: i64 = primes[j]
let v: i64 = i * p
if v > limit { break }
is_comp[v] = 1
if i % p == 0 {
mu[v] = 0
break
}
mu[v] = -mu[i]
j = j + 1
}
}
# Prefix sums of mu (Mertens function)
let pref: ptr<i64> = calloc(limit + 1, 8) as ptr<i64>
let mut s: i64 = 0
for i in 1..(limit + 1) {
s = s + mu[i]
pref[i] = s
}
g_pref = pref
g_limit = limit
g_ht_m_key = calloc(HT_SIZE, 8) as ptr<i64>
g_ht_m_val = calloc(HT_SIZE, 8) as ptr<i64>
g_ht_f_key = calloc(HT_SIZE, 8) as ptr<i64>
g_ht_f_val = calloc(HT_SIZE, 8) as ptr<i64>
let mut ans: i64 = 0
let mut l: i64 = 1
while l <= max_d {
let q: i64 = M / l
let r: i64 = M / q
let mut rr: i64 = r
if rr > max_d { rr = max_d }
let coef: i64 = F_func(rr) - F_func(l - 1)
if coef != 0 {
ans = ans + coef * count_points_nonprimitive(q)
}
l = rr + 1
}
free(mu)
free(is_comp)
free(primes)
free(pref)
free(g_ht_m_key)
free(g_ht_m_val)
free(g_ht_f_key)
free(g_ht_f_val)
return ans
}
function main() -> i32 {
let N: i64 = 1000000000
let M: i64 = 3 * N + 6
let ans: i64 = 2 + 4 * count_points_primitive(M)
printf("%lld\n", ans)
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; }
int64_t floor_sum_i64_i64_i64_i64(int64_t n0, int64_t m0, int64_t a0, int64_t b0);
int64_t icbrt_i64(int64_t n);
int64_t ht_get_ptr_i64_ptr_i64_i64_ptr_i64(int64_t* ht_key, int64_t* ht_val, int64_t key, int64_t* out);
void ht_put_ptr_i64_ptr_i64_i64_i64(int64_t* ht_key, int64_t* ht_val, int64_t key, int64_t val);
int64_t M_func_i64(int64_t n);
int64_t F_func_i64(int64_t n);
int64_t count_points_nonprimitive_i64(int64_t M);
int64_t count_points_primitive_i64(int64_t M);
int32_t main(void);
static const int64_t HT_SIZE = 2097152;
static const int64_t HT_MASK = 2097151;
/* Module statics */
static int64_t* g_pref = NULL;
static int64_t g_limit = 0;
static int64_t* g_ht_m_key = NULL;
static int64_t* g_ht_m_val = NULL;
static int64_t* g_ht_f_key = NULL;
static int64_t* g_ht_f_val = NULL;
int64_t floor_sum_i64_i64_i64_i64(int64_t n0, int64_t m0, int64_t a0, int64_t b0) {
int64_t res = 0;
int64_t n = n0;
int64_t m = m0;
int64_t a = a0;
int64_t b = b0;
while (1) {
if (a >= m) {
res = (res + FLOW_CHECKED_DIV(((((n - 1) * n) * FLOW_CHECKED_DIV((a), (m)))), (2)));
a = FLOW_CHECKED_MOD((a), (m));
}
if (b >= m) {
res = (res + (n * FLOW_CHECKED_DIV((b), (m))));
b = FLOW_CHECKED_MOD((b), (m));
}
int64_t y_max = ((a * n) + b);
if (y_max < m) {
break;
}
b = FLOW_CHECKED_MOD((y_max), (m));
n = FLOW_CHECKED_DIV((y_max), (m));
int64_t tmp = m;
m = a;
a = tmp;
}
return res;
}
int64_t icbrt_i64(int64_t n) {
if (n <= 0) {
return 0;
}
int64_t bits = 0;
int64_t t = n;
while (t != 0) {
bits = (bits + 1);
t = FLOW_CHECKED_SHR((t), (1));
}
int64_t x = FLOW_CHECKED_SHL((1), (FLOW_CHECKED_DIV(((bits + 2)), (3))));
while (1) {
int64_t y = FLOW_CHECKED_DIV((((2 * x) + FLOW_CHECKED_DIV((n), ((x * x))))), (3));
if (y >= x) {
break;
}
x = y;
}
while ((((x + 1) * (x + 1)) * (x + 1)) <= n) {
x = (x + 1);
}
while (((x * x) * x) > n) {
x = (x - 1);
}
return x;
}
int64_t ht_get_ptr_i64_ptr_i64_i64_ptr_i64(int64_t* ht_key, int64_t* ht_val, int64_t key, int64_t* out) {
int64_t h = ((key * 2654435761) & HT_MASK);
while (ht_key[h] != 0) {
if (ht_key[h] == key) {
out[0] = ht_val[h];
return 1;
}
h = ((h + 1) & HT_MASK);
}
return 0;
}
void ht_put_ptr_i64_ptr_i64_i64_i64(int64_t* ht_key, int64_t* ht_val, int64_t key, int64_t val) {
int64_t h = ((key * 2654435761) & HT_MASK);
while (ht_key[h] != 0) {
if (ht_key[h] == key) {
ht_val[h] = val;
return;
}
h = ((h + 1) & HT_MASK);
}
ht_key[h] = key;
ht_val[h] = val;
}
int64_t M_func_i64(int64_t n) {
if (n <= 0) {
return 0;
}
if (n <= g_limit) {
return g_pref[n];
}
int64_t* out = (int64_t*)(((int64_t*)(malloc(8))));
if (ht_get_ptr_i64_ptr_i64_i64_ptr_i64(g_ht_m_key, g_ht_m_val, n, out) == 1) {
int64_t v = out[0];
free(out);
return v;
}
free(out);
int64_t res = 1;
int64_t l = 2;
while (l <= n) {
int64_t q = FLOW_CHECKED_DIV((n), (l));
int64_t r = FLOW_CHECKED_DIV((n), (q));
res = (res - (((r - l) + 1) * M_func_i64(q)));
l = (r + 1);
}
ht_put_ptr_i64_ptr_i64_i64_i64(g_ht_m_key, g_ht_m_val, n, res);
return res;
}
int64_t F_func_i64(int64_t n) {
if (n <= 0) {
return 0;
}
int64_t* out = (int64_t*)(((int64_t*)(malloc(8))));
if (ht_get_ptr_i64_ptr_i64_i64_ptr_i64(g_ht_f_key, g_ht_f_val, n, out) == 1) {
int64_t v = out[0];
free(out);
return v;
}
free(out);
int64_t res = (M_func_i64(n) + F_func_i64(FLOW_CHECKED_DIV((n), (3))));
ht_put_ptr_i64_ptr_i64_i64_i64(g_ht_f_key, g_ht_f_val, n, res);
return res;
}
int64_t count_points_nonprimitive_i64(int64_t M) {
if (M < 28) {
return 0;
}
int64_t n = FLOW_CHECKED_DIV(((M - 18)), (10));
int64_t b = (M - (10 * n));
int64_t total = floor_sum_i64_i64_i64_i64(n, 18, 10, b);
int64_t n3 = FLOW_CHECKED_DIV((n), (3));
int64_t b3 = (M - (30 * n3));
int64_t total3 = floor_sum_i64_i64_i64_i64(n3, 18, 30, b3);
return (total - total3);
}
int64_t count_points_primitive_i64(int64_t M) {
int64_t max_d = FLOW_CHECKED_DIV((M), (28));
if (max_d <= 0) {
return 0;
}
int64_t limit = (icbrt_i64((max_d * max_d)) + 64);
int64_t* mu = (int64_t*)(((int64_t*)(calloc((limit + 1), 8))));
int8_t* is_comp = (int8_t*)(((int8_t*)(calloc((limit + 1), 1))));
int64_t* primes = (int64_t*)(((int64_t*)(malloc(((limit + 1) * 8)))));
int64_t pc = 0;
mu[1] = 1;
int32_t __flow_step_1 = 1;
for (int32_t i = 2; (2 <= (limit + 1)) ? i < (limit + 1) : i > (limit + 1); i += (2 <= (limit + 1)) ? 1 : -1) {
if (is_comp[i] == 0) {
primes[pc] = i;
pc = (pc + 1);
mu[i] = (-1);
}
int64_t j = 0;
while (j < pc) {
int64_t p = primes[j];
int64_t v = (i * p);
if (v > limit) {
break;
}
is_comp[v] = 1;
if (FLOW_CHECKED_MOD((i), (p)) == 0) {
mu[v] = 0;
break;
}
mu[v] = (-mu[i]);
j = (j + 1);
}
}
int64_t* pref = (int64_t*)(((int64_t*)(calloc((limit + 1), 8))));
int64_t s = 0;
int32_t __flow_step_2 = 1;
for (int32_t i = 1; (1 <= (limit + 1)) ? i < (limit + 1) : i > (limit + 1); i += (1 <= (limit + 1)) ? 1 : -1) {
s = (s + mu[i]);
pref[i] = s;
}
g_pref = pref;
g_limit = limit;
g_ht_m_key = ((int64_t*)(calloc(HT_SIZE, 8)));
g_ht_m_val = ((int64_t*)(calloc(HT_SIZE, 8)));
g_ht_f_key = ((int64_t*)(calloc(HT_SIZE, 8)));
g_ht_f_val = ((int64_t*)(calloc(HT_SIZE, 8)));
int64_t ans = 0;
int64_t l = 1;
while (l <= max_d) {
int64_t q = FLOW_CHECKED_DIV((M), (l));
int64_t r = FLOW_CHECKED_DIV((M), (q));
int64_t rr = r;
if (rr > max_d) {
rr = max_d;
}
int64_t coef = (F_func_i64(rr) - F_func_i64((l - 1)));
if (coef != 0) {
ans = (ans + (coef * count_points_nonprimitive_i64(q)));
}
l = (rr + 1);
}
free(mu);
free(is_comp);
free(primes);
free(pref);
free(g_ht_m_key);
free(g_ht_m_val);
free(g_ht_f_key);
free(g_ht_f_val);
return ans;
}
int32_t main(void) {
int64_t N = 1000000000;
int64_t M = ((3 * N) + 6);
int64_t ans = (2 + (4 * count_points_primitive_i64(M)));
printf("%lld\n", ans);
return 0;
}