# Project Euler 915: Giant GCDs
# s(1)=1, s(n+1) = (s(n)-1)^3 + 2
# T(N) = sum_{a=1..N} sum_{b=1..N} gcd(s(s(a)), s(s(b)))
# Compute T(10^8) mod 123456789.
# Uses cycle detection for s(n) mod M, summatory totient via sieve + memoized recursion.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
function malloc(n: i64) -> ptr<void>
}
const MOD: i64 = 123456789
const BASE: i64 = 2000000
const MEMO_CAP: i64 = 1 << 20
let mut s_mod_MOD_arr: ptr<i64> = null
let mut muM: i64 = 0
let mut lamM: i64 = 0
let mut s_mod_lam_arr: ptr<i64> = null
let mut muP: i64 = 0
let mut lamP: i64 = 0
let mut small_exact_arr: ptr<i64> = null
let mut n_small_max: i64 = 0
let mut muM_mod: i64 = 0
let mut phi_prefix_arr: ptr<i64> = null
let mut memo_keys: ptr<i64> = null
let mut memo_vals: ptr<i64> = null
function mulmod(a: i64, b: i64, m: i64) -> i64 {
return ((a as i128) * (b as i128) % (m as i128)) as i64
}
function f_mod(x: i64, m: i64) -> i64 {
let y: i64 = (x + m - 1) % m
let y2: i64 = mulmod(y, y, m)
let y3: i64 = mulmod(y2, y, m)
return (y3 + 2) % m
}
function cycle_info(m: i64) -> ptr<i64> {
let result: ptr<i64> = calloc(2, 8)
let mut tortoise: i64 = f_mod(0, m)
let mut hare: i64 = f_mod(f_mod(0, m), m)
while tortoise != hare {
tortoise = f_mod(tortoise, m)
hare = f_mod(f_mod(hare, m), m)
}
let mut mu: i64 = 0
tortoise = 0
while tortoise != hare {
tortoise = f_mod(tortoise, m)
hare = f_mod(hare, m)
mu = mu + 1
}
let mut lam: i64 = 1
hare = f_mod(tortoise, m)
while tortoise != hare {
hare = f_mod(hare, m)
lam = lam + 1
}
result[0] = mu
result[1] = lam
return result
}
function build_s_mod(m: i64, length: i64) -> ptr<i64> {
let arr: ptr<i64> = calloc(length + 1, 8)
let mut x: i64 = 0
let mut n: i64 = 1
while n <= length {
x = f_mod(x, m)
arr[n] = x
n = n + 1
}
return arr
}
function sieve_phi_prefix(n: i64) -> ptr<i64> {
let phi: ptr<i64> = calloc(n + 1, 8)
let mut i: i64 = 0
while i <= n {
phi[i] = i
i = i + 1
}
i = 2
while i <= n {
if phi[i] == i {
let mut j: i64 = i
while j <= n {
phi[j] = phi[j] - phi[j] / i
j = j + i
}
}
i = i + 1
}
let pref: ptr<i64> = calloc(n + 1, 8)
let mut s: i64 = 0
i = 1
while i <= n {
s = s + phi[i]
pref[i] = s
i = i + 1
}
free(phi as ptr<void>)
return pref
}
function mix64(x: i64) -> i64 {
let mut v: i64 = x
v = v ^ (v >> 30)
v = v * 0xbf58476d1ce4e5b9
v = v ^ (v >> 27)
v = v * 0x94d049bb133111eb
v = v ^ (v >> 31)
return v
}
function s_index_modMOD(k: i64) -> i64 {
if k <= muM + lamM {
return s_mod_MOD_arr[k]
}
let k2: i64 = muM + ((k - muM) % lamM)
return s_mod_MOD_arr[k2]
}
function s_n_mod_lamM(n: i64) -> i64 {
if n <= muP + lamP {
return s_mod_lam_arr[n]
}
let n2: i64 = muP + ((n - muP) % lamP)
return s_mod_lam_arr[n2]
}
function s2_mod(n: i64) -> i64 {
if n <= n_small_max {
return s_index_modMOD(small_exact_arr[n])
}
let k_mod: i64 = s_n_mod_lamM(n)
let diff: i64
if k_mod >= muM_mod {
diff = (k_mod - muM_mod) % lamM
} else {
diff = lamM - ((muM_mod - k_mod) % lamM)
}
if diff == lamM {
return s_mod_MOD_arr[muM]
}
return s_mod_MOD_arr[muM + diff]
}
function phi_sum_memo(n: i64) -> i64 {
if n <= BASE {
return phi_prefix_arr[n]
}
let mut h: i64 = mix64(n) & (MEMO_CAP - 1)
while memo_keys[h] != 0 {
if memo_keys[h] == n { return memo_vals[h] }
h = (h + 1) & (MEMO_CAP - 1)
}
let mut res: i64 = (n * (n + 1)) / 2
let mut l: i64 = 2
while l <= n {
let q: i64 = n / l
let r: i64 = n / q
let sub: i64 = phi_sum_memo(q)
let cnt: i64 = r - l + 1
let val: i128 = (cnt as i128) * (sub as i128)
res = ((res as i128) - val) as i64
l = r + 1
}
memo_keys[h] = n
memo_vals[h] = res
return res
}
function coprime_pairs(m: i64) -> i64 {
let ps: i64 = phi_sum_memo(m) % MOD
return (2 * ps - 1 + MOD) % MOD
}
function prefix_s2(n: i64, start: i64, period: i64, small_prefix: ptr<i64>, period_prefix: ptr<i64>, period_sum: i64) -> i64 {
if n == 0 { return 0 }
if n < start { return small_prefix[n] }
let base: i64 = 0
if start > 0 { base = small_prefix[start - 1] }
let t: i64 = n - (start - 1)
let full: i64 = t / period
let rem: i64 = t % period
return (base + mulmod(full % MOD, period_sum, MOD) + period_prefix[rem]) % MOD
}
function main() -> i32 {
let N: i64 = 100000000
# Periodicity of s(n) mod MOD
let ci: ptr<i64> = cycle_info(MOD)
muM = ci[0]
lamM = ci[1]
free(ci as ptr<void>)
s_mod_MOD_arr = build_s_mod(MOD, muM + lamM)
# Periodicity of s(n) mod lamM
let ci2: ptr<i64> = cycle_info(lamM)
muP = ci2[0]
lamP = ci2[1]
free(ci2 as ptr<void>)
s_mod_lam_arr = build_s_mod(lamM, muP + lamP)
# Compute exact s(n) for small n until s(n) > muM
small_exact_arr = calloc(20, 8)
small_exact_arr[0] = 0
let mut x: i64 = 0
let mut n: i64 = 0
while true {
n = n + 1
let xm1: i128 = (x as i128) - 1
let xc: i128 = xm1 * xm1 * xm1 + 2
x = xc as i64
small_exact_arr[n] = x
if x > muM { break }
}
n_small_max = n - 1
muM_mod = muM % lamM
# s2_mod becomes periodic once s(n) mod lamM is in its cycle
let mut start: i64 = muP
if n_small_max + 1 > start { start = n_small_max + 1 }
if 1 > start { start = 1 }
let period: i64 = lamP
let period_vals: ptr<i64> = calloc(period, 8)
let mut i: i64 = 0
while i < period {
period_vals[i] = s2_mod(start + i) % MOD
i = i + 1
}
let small_prefix: ptr<i64> = calloc(start, 8)
let mut acc: i64 = 0
i = 1
while i < start {
acc = (acc + s2_mod(i)) % MOD
small_prefix[i] = acc
i = i + 1
}
let period_prefix: ptr<i64> = calloc(period + 1, 8)
let mut accp: i64 = 0
i = 0
while i < period {
accp = (accp + period_vals[i]) % MOD
period_prefix[i + 1] = accp
i = i + 1
}
let period_sum: i64 = period_prefix[period]
# Summatory totient
phi_prefix_arr = sieve_phi_prefix(BASE)
memo_keys = calloc(MEMO_CAP, 8)
memo_vals = calloc(MEMO_CAP, 8)
# Block over d where floor(N/d) is constant
let mut ans: i64 = 0
let mut l: i64 = 1
while l <= N {
let q: i64 = N / l
let r: i64 = N / q
let pr: i64 = prefix_s2(r, start, period, small_prefix, period_prefix, period_sum)
let pl: i64 = prefix_s2(l - 1, start, period, small_prefix, period_prefix, period_sum)
let sum_s2: i64 = (pr + MOD - pl) % MOD
ans = (ans + mulmod(sum_s2, coprime_pairs(q), MOD)) % MOD
l = r + 1
}
printf("%lld\n", ans % MOD)
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 mulmod_i64_i64_i64(int64_t a, int64_t b, int64_t m);
int64_t f_mod_i64_i64(int64_t x, int64_t m);
int64_t* cycle_info_i64(int64_t m);
int64_t* build_s_mod_i64_i64(int64_t m, int64_t length);
int64_t* sieve_phi_prefix_i64(int64_t n);
int64_t mix64_i64(int64_t x);
int64_t s_index_modMOD_i64(int64_t k);
int64_t s_n_mod_lamM_i64(int64_t n);
int64_t s2_mod_i64(int64_t n);
int64_t phi_sum_memo_i64(int64_t n);
int64_t coprime_pairs_i64(int64_t m);
int64_t prefix_s2_i64_i64_i64_ptr_i64_ptr_i64_i64(int64_t n, int64_t start, int64_t period, int64_t* small_prefix, int64_t* period_prefix, int64_t period_sum);
int32_t main(void);
static const int64_t MOD = 123456789;
static const int64_t BASE = 2000000;
static const int64_t MEMO_CAP = FLOW_CHECKED_SHL((1), (20));
/* Module statics */
static int64_t* s_mod_MOD_arr = NULL;
static int64_t muM = 0;
static int64_t lamM = 0;
static int64_t* s_mod_lam_arr = NULL;
static int64_t muP = 0;
static int64_t lamP = 0;
static int64_t* small_exact_arr = NULL;
static int64_t n_small_max = 0;
static int64_t muM_mod = 0;
static int64_t* phi_prefix_arr = NULL;
static int64_t* memo_keys = NULL;
static int64_t* memo_vals = NULL;
int64_t mulmod_i64_i64_i64(int64_t a, int64_t b, int64_t m) {
return ((int64_t)(FLOW_CHECKED_MOD(((((__int128)(a)) * ((__int128)(b)))), (((__int128)(m))))));
}
int64_t f_mod_i64_i64(int64_t x, int64_t m) {
int64_t y = FLOW_CHECKED_MOD((((x + m) - 1)), (m));
int64_t y2 = mulmod_i64_i64_i64(y, y, m);
int64_t y3 = mulmod_i64_i64_i64(y2, y, m);
return FLOW_CHECKED_MOD(((y3 + 2)), (m));
}
int64_t* cycle_info_i64(int64_t m) {
int64_t* result = (int64_t*)(calloc(2, 8));
int64_t tortoise = f_mod_i64_i64(0, m);
int64_t hare = f_mod_i64_i64(f_mod_i64_i64(0, m), m);
while (tortoise != hare) {
tortoise = f_mod_i64_i64(tortoise, m);
hare = f_mod_i64_i64(f_mod_i64_i64(hare, m), m);
}
int64_t mu = 0;
tortoise = 0;
while (tortoise != hare) {
tortoise = f_mod_i64_i64(tortoise, m);
hare = f_mod_i64_i64(hare, m);
mu = (mu + 1);
}
int64_t lam = 1;
hare = f_mod_i64_i64(tortoise, m);
while (tortoise != hare) {
hare = f_mod_i64_i64(hare, m);
lam = (lam + 1);
}
result[0] = mu;
result[1] = lam;
return result;
}
int64_t* build_s_mod_i64_i64(int64_t m, int64_t length) {
int64_t* arr = (int64_t*)(calloc((length + 1), 8));
int64_t x = 0;
int64_t n = 1;
while (n <= length) {
x = f_mod_i64_i64(x, m);
arr[n] = x;
n = (n + 1);
}
return arr;
}
int64_t* sieve_phi_prefix_i64(int64_t n) {
int64_t* phi = (int64_t*)(calloc((n + 1), 8));
int64_t i = 0;
while (i <= n) {
phi[i] = i;
i = (i + 1);
}
i = 2;
while (i <= n) {
if (phi[i] == i) {
int64_t j = i;
while (j <= n) {
phi[j] = (phi[j] - FLOW_CHECKED_DIV((phi[j]), (i)));
j = (j + i);
}
}
i = (i + 1);
}
int64_t* pref = (int64_t*)(calloc((n + 1), 8));
int64_t s = 0;
i = 1;
while (i <= n) {
s = (s + phi[i]);
pref[i] = s;
i = (i + 1);
}
free(((void*)(phi)));
return pref;
}
int64_t mix64_i64(int64_t x) {
int64_t v = x;
v = (v ^ FLOW_CHECKED_SHR((v), (30)));
v = (v * ((__int128)0xBF58476D1CE4E5B9ULL));
v = (v ^ FLOW_CHECKED_SHR((v), (27)));
v = (v * ((__int128)0x94D049BB133111EBULL));
v = (v ^ FLOW_CHECKED_SHR((v), (31)));
return v;
}
int64_t s_index_modMOD_i64(int64_t k) {
if (k <= (muM + lamM)) {
return s_mod_MOD_arr[k];
}
int64_t k2 = (muM + FLOW_CHECKED_MOD(((k - muM)), (lamM)));
return s_mod_MOD_arr[k2];
}
int64_t s_n_mod_lamM_i64(int64_t n) {
if (n <= (muP + lamP)) {
return s_mod_lam_arr[n];
}
int64_t n2 = (muP + FLOW_CHECKED_MOD(((n - muP)), (lamP)));
return s_mod_lam_arr[n2];
}
int64_t s2_mod_i64(int64_t n) {
if (n <= n_small_max) {
return s_index_modMOD_i64(small_exact_arr[n]);
}
int64_t k_mod = s_n_mod_lamM_i64(n);
int64_t diff;
if (k_mod >= muM_mod) {
diff = FLOW_CHECKED_MOD(((k_mod - muM_mod)), (lamM));
} else {
diff = (lamM - FLOW_CHECKED_MOD(((muM_mod - k_mod)), (lamM)));
}
if (diff == lamM) {
return s_mod_MOD_arr[muM];
}
return s_mod_MOD_arr[(muM + diff)];
}
int64_t phi_sum_memo_i64(int64_t n) {
if (n <= BASE) {
return phi_prefix_arr[n];
}
int64_t h = (mix64_i64(n) & (MEMO_CAP - 1));
while (memo_keys[h] != 0) {
if (memo_keys[h] == n) {
return memo_vals[h];
}
h = ((h + 1) & (MEMO_CAP - 1));
}
int64_t res = FLOW_CHECKED_DIV(((n * (n + 1))), (2));
int64_t l = 2;
while (l <= n) {
int64_t q = FLOW_CHECKED_DIV((n), (l));
int64_t r = FLOW_CHECKED_DIV((n), (q));
int64_t sub = phi_sum_memo_i64(q);
int64_t cnt = ((r - l) + 1);
__int128 val = (((__int128)(cnt)) * ((__int128)(sub)));
res = ((int64_t)((((__int128)(res)) - val)));
l = (r + 1);
}
memo_keys[h] = n;
memo_vals[h] = res;
return res;
}
int64_t coprime_pairs_i64(int64_t m) {
int64_t ps = FLOW_CHECKED_MOD((phi_sum_memo_i64(m)), (MOD));
return FLOW_CHECKED_MOD(((((2 * ps) - 1) + MOD)), (MOD));
}
int64_t prefix_s2_i64_i64_i64_ptr_i64_ptr_i64_i64(int64_t n, int64_t start, int64_t period, int64_t* small_prefix, int64_t* period_prefix, int64_t period_sum) {
if (n == 0) {
return 0;
}
if (n < start) {
return small_prefix[n];
}
int64_t base = 0;
if (start > 0) {
base = small_prefix[(start - 1)];
}
int64_t t = (n - (start - 1));
int64_t full = FLOW_CHECKED_DIV((t), (period));
int64_t rem = FLOW_CHECKED_MOD((t), (period));
return FLOW_CHECKED_MOD((((base + mulmod_i64_i64_i64(FLOW_CHECKED_MOD((full), (MOD)), period_sum, MOD)) + period_prefix[rem])), (MOD));
}
int32_t main(void) {
int64_t N = 100000000;
int64_t* ci = (int64_t*)(cycle_info_i64(MOD));
muM = ci[0];
lamM = ci[1];
free(((void*)(ci)));
s_mod_MOD_arr = build_s_mod_i64_i64(MOD, (muM + lamM));
int64_t* ci2 = (int64_t*)(cycle_info_i64(lamM));
muP = ci2[0];
lamP = ci2[1];
free(((void*)(ci2)));
s_mod_lam_arr = build_s_mod_i64_i64(lamM, (muP + lamP));
small_exact_arr = calloc(20, 8);
small_exact_arr[0] = 0;
int64_t x = 0;
int64_t n = 0;
while (1) {
n = (n + 1);
__int128 xm1 = (((__int128)(x)) - 1);
__int128 xc = (((xm1 * xm1) * xm1) + 2);
x = ((int64_t)(xc));
small_exact_arr[n] = x;
if (x > muM) {
break;
}
}
n_small_max = (n - 1);
muM_mod = FLOW_CHECKED_MOD((muM), (lamM));
int64_t start = muP;
if ((n_small_max + 1) > start) {
start = (n_small_max + 1);
}
if (1 > start) {
start = 1;
}
int64_t period = lamP;
int64_t* period_vals = (int64_t*)(calloc(period, 8));
int64_t i = 0;
while (i < period) {
period_vals[i] = FLOW_CHECKED_MOD((s2_mod_i64((start + i))), (MOD));
i = (i + 1);
}
int64_t* small_prefix = (int64_t*)(calloc(start, 8));
int64_t acc = 0;
i = 1;
while (i < start) {
acc = FLOW_CHECKED_MOD(((acc + s2_mod_i64(i))), (MOD));
small_prefix[i] = acc;
i = (i + 1);
}
int64_t* period_prefix = (int64_t*)(calloc((period + 1), 8));
int64_t accp = 0;
i = 0;
while (i < period) {
accp = FLOW_CHECKED_MOD(((accp + period_vals[i])), (MOD));
period_prefix[(i + 1)] = accp;
i = (i + 1);
}
int64_t period_sum = period_prefix[period];
phi_prefix_arr = sieve_phi_prefix_i64(BASE);
memo_keys = calloc(MEMO_CAP, 8);
memo_vals = calloc(MEMO_CAP, 8);
int64_t ans = 0;
int64_t l = 1;
while (l <= N) {
int64_t q = FLOW_CHECKED_DIV((N), (l));
int64_t r = FLOW_CHECKED_DIV((N), (q));
int64_t pr = prefix_s2_i64_i64_i64_ptr_i64_ptr_i64_i64(r, start, period, small_prefix, period_prefix, period_sum);
int64_t pl = prefix_s2_i64_i64_i64_ptr_i64_ptr_i64_i64((l - 1), start, period, small_prefix, period_prefix, period_sum);
int64_t sum_s2 = FLOW_CHECKED_MOD((((pr + MOD) - pl)), (MOD));
ans = FLOW_CHECKED_MOD(((ans + mulmod_i64_i64_i64(sum_s2, coprime_pairs_i64(q), MOD))), (MOD));
l = (r + 1);
}
printf("%lld\n", FLOW_CHECKED_MOD((ans), (MOD)));
return 0;
}