# Project Euler 241
# Sum of n <= 10^18 with sigma(n)/n a half-integer.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
function gcd(a0: i64, b0: i64) -> i64 {
let mut a: i64 = a0
let mut b: i64 = b0
while b != 0 {
let t: i64 = a % b
a = b
b = t
}
return a
}
function addmod(a: i64, b: i64, mod: i64) -> i64 {
let s: i64 = a + b
if s >= mod || s < a { return s - mod }
return s
}
function mulmod(a0: i64, b0: i64, mod: i64) -> i64 {
let mut a: i64 = a0 % mod
let mut b: i64 = b0 % mod
if a < 0 { a = a + mod }
if b < 0 { b = b + mod }
let mut result: i64 = 0
while b > 0 {
if b % 2 == 1 { result = addmod(result, a, mod) }
a = addmod(a, a, mod)
b = b / 2
}
return result
}
function modpow(base: i64, exp: i64, mod: i64) -> i64 {
let mut r: i64 = 1
let mut b: i64 = base % mod
let mut e: i64 = exp
while e > 0 {
if e % 2 == 1 { r = mulmod(r, b, mod) }
b = mulmod(b, b, mod)
e = e / 2
}
return r
}
function is_prime(value: i64) -> bool {
if value < 2 { return false }
if value % 2 == 0 { return value == 2 }
if value % 3 == 0 { return value == 3 }
if value % 5 == 0 { return value == 5 }
if value % 7 == 0 { return value == 7 }
if value % 11 == 0 { return value == 11 }
if value % 13 == 0 { return value == 13 }
if value % 17 == 0 { return value == 17 }
if value % 19 == 0 { return value == 19 }
if value % 23 == 0 { return value == 23 }
if value % 29 == 0 { return value == 29 }
let mut d: i64 = value - 1
let mut shifts: i32 = 0
while d % 2 == 0 { d = d / 2; shifts = shifts + 1 }
let bases: ptr<i64> = calloc(7, 8)
bases[0] = 2; bases[1] = 325; bases[2] = 9375; bases[3] = 28178
bases[4] = 450775; bases[5] = 9780504; bases[6] = 1795265022
for bi in 0..7 {
let base: i64 = bases[bi] % value
if base != 0 {
let mut x: i64 = modpow(base, d, value)
if !(x == 1 || x == value - 1) {
let mut ok: bool = false
let mut r: i32 = 0
while r < shifts - 1 {
x = mulmod(x, x, value)
if x == value - 1 { ok = true; break }
r = r + 1
}
if !ok { free(bases); return false }
}
}
}
free(bases)
return true
}
function pollard_rho(n: i64) -> i64 {
if n % 2 == 0 { return 2 }
let mut c: i64 = 1
while true {
let mut x: i64 = 2
let mut y: i64 = 2
let mut d: i64 = 1
while d == 1 {
x = addmod(mulmod(x, x, n), c, n)
y = addmod(mulmod(y, y, n), c, n)
y = addmod(mulmod(y, y, n), c, n)
let mut diff: i64 = x - y
if diff < 0 { diff = 0 - diff }
d = gcd(diff, n)
}
if d != n { return d }
c = c + 1
}
return n
}
# Factor into primes array (with multiplicity), then group
function factor_smallest(n0: i64) -> i64 {
# return smallest prime factor
let mut n: i64 = n0
if n % 2 == 0 { return 2 }
if is_prime(n) { return n }
let f: i64 = pollard_rho(n)
while !is_prime(f) {
f = pollard_rho(f)
}
# also check n/f factors - we need smallest prime factor of n0
let mut spf: i64 = f
let mut rem: i64 = n0
# fully factor via trial + pollard
let mut p: i64 = 2
while p * p <= rem && p <= 1000000 {
if rem % p == 0 { return p }
if p == 2 { p = 3 } else { p = p + 2 }
}
if rem > 1 && is_prime(rem) { return rem }
# pollard tree
let mut stack: ptr<i64> = calloc(64, 8)
let mut sp: i32 = 0
stack[0] = n0; sp = 1
let mut best: i64 = n0
while sp > 0 {
sp = sp - 1
let cur: i64 = stack[sp]
if cur < best && is_prime(cur) { best = cur; continue }
if is_prime(cur) {
if cur < best { best = cur }
continue
}
let g: i64 = pollard_rho(cur)
stack[sp] = g; sp = sp + 1
stack[sp] = cur / g; sp = sp + 1
}
free(stack)
return best
}
let mut TOTAL: i64 = 0
const LIMIT: i64 = 1000000000000000000
function search_target(target_num: i64) -> void {
# stack of (n, numerator, denominator) - ignore used_primes by checking factor of denom
let CAP: i64 = 100000
let sn: ptr<i64> = calloc(CAP, 8)
let snum: ptr<i64> = calloc(CAP, 8)
let sden: ptr<i64> = calloc(CAP, 8)
let mut sp: i64 = 0
sn[0] = 1; snum[0] = target_num; sden[0] = 2; sp = 1
while sp > 0 {
sp = sp - 1
let n: i64 = sn[sp]
let numerator: i64 = snum[sp]
let denominator: i64 = sden[sp]
if numerator == denominator {
TOTAL = TOTAL + n
continue
}
if numerator < denominator { continue }
if n > LIMIT / denominator { continue }
if denominator == 1 { continue }
let p: i64 = factor_smallest(denominator)
# min exponent: how many times p divides denominator? for the forced min, use valuation of denom
# From python: p, min_exponent = factor_items_tuple(denominator)[0]
# factor_items_tuple returns sorted unique primes with exponents of the full factorization
# So we need the smallest prime and its exponent in denominator
let mut min_exp: i32 = 0
let mut tmp: i64 = denominator
while tmp % p == 0 {
tmp = tmp / p
min_exp = min_exp + 1
}
# check p not already used: if n % p == 0 then used
if n % p == 0 { continue }
let mut prime_power: i64 = 1
let mut sigma_power: i64 = 1
for e in 0..min_exp {
prime_power = prime_power * p
sigma_power = sigma_power + prime_power
}
while n <= LIMIT / prime_power {
let mut new_num: i64 = numerator * prime_power
let mut new_den: i64 = denominator * sigma_power
let g: i64 = gcd(new_num, new_den)
new_num = new_num / g
new_den = new_den / g
if new_num < new_den { break }
sn[sp] = n * prime_power
snum[sp] = new_num
sden[sp] = new_den
sp = sp + 1
if prime_power > LIMIT / p { break }
prime_power = prime_power * p
sigma_power = sigma_power + prime_power
}
}
free(sn); free(snum); free(sden)
}
function main() -> i32 {
let targets: ptr<i64> = calloc(6, 8)
targets[0] = 3; targets[1] = 5; targets[2] = 7
targets[3] = 9; targets[4] = 11; targets[5] = 13
for i in 0..6 {
search_target(targets[i])
}
printf("%lld\n", TOTAL)
free(targets)
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 gcd_i64_i64(int64_t a0, int64_t b0);
int64_t addmod_i64_i64_i64(int64_t a, int64_t b, int64_t mod);
int64_t mulmod_i64_i64_i64(int64_t a0, int64_t b0, int64_t mod);
int64_t modpow_i64_i64_i64(int64_t base, int64_t exp, int64_t mod);
bool is_prime_i64(int64_t value);
int64_t pollard_rho_i64(int64_t n);
int64_t factor_smallest_i64(int64_t n0);
void search_target_i64(int64_t target_num);
int32_t main(void);
static const int64_t LIMIT = 1000000000000000000;
/* Module statics */
static int64_t TOTAL = 0;
int64_t gcd_i64_i64(int64_t a0, int64_t b0) {
int64_t a = a0;
int64_t b = b0;
while (b != 0) {
int64_t t = FLOW_CHECKED_MOD((a), (b));
a = b;
b = t;
}
return a;
}
int64_t addmod_i64_i64_i64(int64_t a, int64_t b, int64_t mod) {
int64_t s = (a + b);
if ((s >= mod || s < a)) {
return (s - mod);
}
return s;
}
int64_t mulmod_i64_i64_i64(int64_t a0, int64_t b0, int64_t mod) {
int64_t a = FLOW_CHECKED_MOD((a0), (mod));
int64_t b = FLOW_CHECKED_MOD((b0), (mod));
if (a < 0) {
a = (a + mod);
}
if (b < 0) {
b = (b + mod);
}
int64_t result = 0;
while (b > 0) {
if (FLOW_CHECKED_MOD((b), (2)) == 1) {
result = addmod_i64_i64_i64(result, a, mod);
}
a = addmod_i64_i64_i64(a, a, mod);
b = FLOW_CHECKED_DIV((b), (2));
}
return result;
}
int64_t modpow_i64_i64_i64(int64_t base, int64_t exp, int64_t mod) {
int64_t r = 1;
int64_t b = FLOW_CHECKED_MOD((base), (mod));
int64_t e = exp;
while (e > 0) {
if (FLOW_CHECKED_MOD((e), (2)) == 1) {
r = mulmod_i64_i64_i64(r, b, mod);
}
b = mulmod_i64_i64_i64(b, b, mod);
e = FLOW_CHECKED_DIV((e), (2));
}
return r;
}
bool is_prime_i64(int64_t value) {
if (value < 2) {
return 0;
}
if (FLOW_CHECKED_MOD((value), (2)) == 0) {
return value == 2;
}
if (FLOW_CHECKED_MOD((value), (3)) == 0) {
return value == 3;
}
if (FLOW_CHECKED_MOD((value), (5)) == 0) {
return value == 5;
}
if (FLOW_CHECKED_MOD((value), (7)) == 0) {
return value == 7;
}
if (FLOW_CHECKED_MOD((value), (11)) == 0) {
return value == 11;
}
if (FLOW_CHECKED_MOD((value), (13)) == 0) {
return value == 13;
}
if (FLOW_CHECKED_MOD((value), (17)) == 0) {
return value == 17;
}
if (FLOW_CHECKED_MOD((value), (19)) == 0) {
return value == 19;
}
if (FLOW_CHECKED_MOD((value), (23)) == 0) {
return value == 23;
}
if (FLOW_CHECKED_MOD((value), (29)) == 0) {
return value == 29;
}
int64_t d = (value - 1);
int32_t shifts = 0;
while (FLOW_CHECKED_MOD((d), (2)) == 0) {
d = FLOW_CHECKED_DIV((d), (2));
shifts = (shifts + 1);
}
int64_t* bases = (int64_t*)(calloc(7, 8));
bases[0] = 2;
bases[1] = 325;
bases[2] = 9375;
bases[3] = 28178;
bases[4] = 450775;
bases[5] = 9780504;
bases[6] = 1795265022;
int32_t __flow_step_1 = 1;
for (int32_t bi = 0; (0 <= 7) ? bi < 7 : bi > 7; bi += (0 <= 7) ? 1 : -1) {
int64_t base = FLOW_CHECKED_MOD((bases[bi]), (value));
if (base != 0) {
int64_t x = modpow_i64_i64_i64(base, d, value);
if ((!((x == 1 || x == (value - 1))))) {
bool ok = 0;
int32_t r = 0;
while (r < (shifts - 1)) {
x = mulmod_i64_i64_i64(x, x, value);
if (x == (value - 1)) {
ok = 1;
break;
}
r = (r + 1);
}
if ((!(ok))) {
free(bases);
return 0;
}
}
}
}
free(bases);
return 1;
}
int64_t pollard_rho_i64(int64_t n) {
if (FLOW_CHECKED_MOD((n), (2)) == 0) {
return 2;
}
int64_t c = 1;
while (1) {
int64_t x = 2;
int64_t y = 2;
int64_t d = 1;
while (d == 1) {
x = addmod_i64_i64_i64(mulmod_i64_i64_i64(x, x, n), c, n);
y = addmod_i64_i64_i64(mulmod_i64_i64_i64(y, y, n), c, n);
y = addmod_i64_i64_i64(mulmod_i64_i64_i64(y, y, n), c, n);
int64_t diff = (x - y);
if (diff < 0) {
diff = (0 - diff);
}
d = gcd_i64_i64(diff, n);
}
if (d != n) {
return d;
}
c = (c + 1);
}
return n;
}
int64_t factor_smallest_i64(int64_t n0) {
int64_t n = n0;
if (FLOW_CHECKED_MOD((n), (2)) == 0) {
return 2;
}
if (is_prime_i64(n)) {
return n;
}
int64_t f = pollard_rho_i64(n);
while ((!(is_prime_i64(f)))) {
f = pollard_rho_i64(f);
}
int64_t spf = f;
int64_t rem = n0;
int64_t p = 2;
while (((p * p) <= rem && p <= 1000000)) {
if (FLOW_CHECKED_MOD((rem), (p)) == 0) {
return p;
}
if (p == 2) {
p = 3;
} else {
p = (p + 2);
}
}
if ((rem > 1 && is_prime_i64(rem))) {
return rem;
}
int64_t* stack = (int64_t*)(calloc(64, 8));
int32_t sp = 0;
stack[0] = n0;
sp = 1;
int64_t best = n0;
while (sp > 0) {
sp = (sp - 1);
int64_t cur = stack[sp];
if ((cur < best && is_prime_i64(cur))) {
best = cur;
continue;
}
if (is_prime_i64(cur)) {
if (cur < best) {
best = cur;
}
continue;
}
int64_t g = pollard_rho_i64(cur);
stack[sp] = g;
sp = (sp + 1);
stack[sp] = FLOW_CHECKED_DIV((cur), (g));
sp = (sp + 1);
}
free(stack);
return best;
}
void search_target_i64(int64_t target_num) {
int64_t CAP = 100000;
int64_t* sn = (int64_t*)(calloc(CAP, 8));
int64_t* snum = (int64_t*)(calloc(CAP, 8));
int64_t* sden = (int64_t*)(calloc(CAP, 8));
int64_t sp = 0;
sn[0] = 1;
snum[0] = target_num;
sden[0] = 2;
sp = 1;
while (sp > 0) {
sp = (sp - 1);
int64_t n = sn[sp];
int64_t numerator = snum[sp];
int64_t denominator = sden[sp];
if (numerator == denominator) {
TOTAL = (TOTAL + n);
continue;
}
if (numerator < denominator) {
continue;
}
if (n > FLOW_CHECKED_DIV((LIMIT), (denominator))) {
continue;
}
if (denominator == 1) {
continue;
}
int64_t p = factor_smallest_i64(denominator);
int32_t min_exp = 0;
int64_t tmp = denominator;
while (FLOW_CHECKED_MOD((tmp), (p)) == 0) {
tmp = FLOW_CHECKED_DIV((tmp), (p));
min_exp = (min_exp + 1);
}
if (FLOW_CHECKED_MOD((n), (p)) == 0) {
continue;
}
int64_t prime_power = 1;
int64_t sigma_power = 1;
int32_t __flow_step_2 = 1;
for (int32_t e = 0; (0 <= min_exp) ? e < min_exp : e > min_exp; e += (0 <= min_exp) ? 1 : -1) {
prime_power = (prime_power * p);
sigma_power = (sigma_power + prime_power);
}
while (n <= FLOW_CHECKED_DIV((LIMIT), (prime_power))) {
int64_t new_num = (numerator * prime_power);
int64_t new_den = (denominator * sigma_power);
int64_t g = gcd_i64_i64(new_num, new_den);
new_num = FLOW_CHECKED_DIV((new_num), (g));
new_den = FLOW_CHECKED_DIV((new_den), (g));
if (new_num < new_den) {
break;
}
sn[sp] = (n * prime_power);
snum[sp] = new_num;
sden[sp] = new_den;
sp = (sp + 1);
if (prime_power > FLOW_CHECKED_DIV((LIMIT), (p))) {
break;
}
prime_power = (prime_power * p);
sigma_power = (sigma_power + prime_power);
}
}
free(sn);
free(snum);
free(sden);
}
int32_t main(void) {
int64_t* targets = (int64_t*)(calloc(6, 8));
targets[0] = 3;
targets[1] = 5;
targets[2] = 7;
targets[3] = 9;
targets[4] = 11;
targets[5] = 13;
int32_t __flow_step_3 = 1;
for (int32_t i = 0; (0 <= 6) ? i < 6 : i > 6; i += (0 <= 6) ? 1 : -1) {
search_target_i64(targets[i]);
}
printf("%lld\n", TOTAL);
free(targets);
return 0;
}