# Project Euler 801: x^y = y^x (mod p).
# Sum of f(p) for primes p in [10^16, 10^16 + 10^6], mod 993353399.
# Pure Flow port of the native C solver. Uses i128 for mulmod/powmod.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
const MOD: i64 = 993353399
# 64-bit modular multiply via i128.
function mulmod(a: i64, b: i64, m: i64) -> i64 {
return ((a as i128) * (b as i128) % (m as i128)) as i64
}
function powmod(a0: i64, e0: i64, m: i64) -> i64 {
let mut r: i64 = 1 % m
let mut a: i64 = a0 % m
let mut e: i64 = e0
while e != 0 {
if e % 2 == 1 {
r = mulmod(r, a, m)
}
a = mulmod(a, a, m)
e = e / 2
}
return r
}
# Deterministic Miller-Rabin for 64-bit.
function is_prime_64(n: i64) -> i32 {
if n < 2 {
return 0
}
# small primes trial
let sp: array<i64, 12> = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
let mut i: i64 = 0
while i < 12 {
if n == sp[i] {
return 1
}
if n % sp[i] == 0 {
return 0
}
i = i + 1
}
let mut d: i64 = n - 1
let mut s: i64 = 0
while d % 2 == 0 {
d = d / 2
s = s + 1
}
let bases: array<i64, 7> = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
let mut bi: i64 = 0
while bi < 7 {
let a: i64 = bases[bi]
if a % n == 0 {
bi = bi + 1
} else {
let mut x: i64 = powmod(a, d, n)
if x == 1 || x == n - 1 {
bi = bi + 1
} else {
let mut comp: i32 = 1
let mut j: i64 = 0
while j < s - 1 {
x = mulmod(x, x, n)
if x == n - 1 {
comp = 0
j = s
}
j = j + 1
}
if comp != 0 {
return 0
}
bi = bi + 1
}
}
}
return 1
}
function gcd_u64(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
}
let mut rng_state: i64 = -7046029254386353131
function rand64() -> i64 {
let mut x: i64 = rng_state
x = x ^ (x >> 12)
x = x ^ (x << 25)
x = x ^ (x >> 27)
rng_state = x
return x * 2685821657736338717
}
function pollard_rho(n: i64) -> i64 {
if n % 2 == 0 {
return 2
}
if n % 3 == 0 {
return 3
}
while true {
let c: i64 = rand64() % (n - 1) + 1
let mut x: i64 = rand64() % (n - 2) + 2
let mut y: i64 = x
let mut d: i64 = 1
while d == 1 {
x = (mulmod(x, x, n) + c) % n
y = (mulmod(y, y, n) + c) % n
y = (mulmod(y, y, n) + c) % n
let mut diff: i64 = x - y
if diff < 0 {
diff = -diff
}
d = gcd_u64(diff, n)
}
if d != n {
return d
}
}
return 0
}
# Factorize n into prime/exp arrays. Returns count.
function factorize(n0: i64, fp: ptr<i64>, ep: ptr<i64>) -> i64 {
let mut outn: i64 = 0
let mut n: i64 = n0
if n <= 1 {
return 0
}
let mut p: i64 = 2
while p * p <= n && p <= 10000 {
if n % p == 0 {
let mut e: i64 = 0
while n % p == 0 {
n = n / p
e = e + 1
}
fp[outn] = p
ep[outn] = e
outn = outn + 1
}
p = p + 1
}
if n == 1 {
return outn
}
if is_prime_64(n) != 0 {
fp[outn] = n
ep[outn] = 1
outn = outn + 1
return outn
}
# Pollard-Rho with explicit stack
let stack: ptr<i64> = calloc(64, 8) as ptr<i64>
let mut sp: i64 = 0
stack[sp] = n
sp = sp + 1
while sp > 0 {
sp = sp - 1
let m: i64 = stack[sp]
if m == 1 {
} else {
if is_prime_64(m) != 0 {
let mut found: i32 = 0
let mut i: i64 = 0
while i < outn {
if fp[i] == m {
ep[i] = ep[i] + 1
found = 1
i = outn
}
i = i + 1
}
if found == 0 {
fp[outn] = m
ep[outn] = 1
outn = outn + 1
}
} else {
let d: i64 = pollard_rho(m)
stack[sp] = d
sp = sp + 1
stack[sp] = m / d
sp = sp + 1
}
}
}
free(stack)
return outn
}
# powmod for small modulus (fits in 32 bits).
function powmod_small(a0: i64, e0: i64, m: i64) -> i64 {
let mut r: i64 = 1 % m
let mut a: i64 = a0 % m
if a < 0 {
a = a + m
}
let mut e: i64 = e0
while e > 0 {
if e % 2 == 1 {
r = r * a % m
}
a = a * a % m
e = e / 2
}
return r
}
# g(q^e) mod.
function g_prime_power_mod(q: i64, e: i64, modv: i64) -> i64 {
let qm: i64 = q % modv
let mut s: i64 = 0
let mut t: i64 = 1
while t <= e {
let expv: i64 = 3 * e - t - 2
let term: i64 = powmod_small(qm, expv, modv)
s = (s + (t * t) % modv * term) % modv
t = t + 1
}
let qm1: i64 = (q - 1) % modv
let term1: i64 = powmod_small(qm1, 3, modv) * s % modv
let inner: i64 = (e * (q - 1) + q) % modv
let term2: i64 = powmod_small(qm, 2 * e - 2, modv) * powmod_small(inner, 2, modv) % modv
return (term1 + term2) % modv
}
function g_from_factorization(fp: ptr<i64>, ep: ptr<i64>, nf: i64, modv: i64) -> i64 {
let mut g: i64 = 1 % modv
let mut i: i64 = 0
while i < nf {
g = g * g_prime_power_mod(fp[i], ep[i], modv) % modv
i = i + 1
}
return g
}
function f_of_prime(p: i64, modv: i64) -> i64 {
let m: i64 = p - 1
let fp: ptr<i64> = calloc(64, 8) as ptr<i64>
let ep: ptr<i64> = calloc(64, 8) as ptr<i64>
let nf: i64 = factorize(m, fp, ep)
let g: i64 = g_from_factorization(fp, ep, nf, modv)
free(fp)
free(ep)
let mm: i64 = m % modv
return (mm * mm + g) % modv
}
function main() -> i32 {
let modv: i64 = MOD
let lo: i64 = 10000000000000000
let hi: i64 = lo + 1000000
let length: i64 = hi - lo + 1
let is_comp: ptr<i8> = calloc(length, 1) as ptr<i8>
# sieve small primes up to 200000
let limit: i64 = 200000
let sieve: ptr<i8> = calloc(limit + 1, 1) as ptr<i8>
let primes: ptr<i64> = calloc(20000, 8) as ptr<i64>
let mut pc: i64 = 0
let mut pp: i64 = 2
while pp <= limit {
if sieve[pp] == 0 {
primes[pc] = pp
pc = pc + 1
let mut j: i64 = pp * pp
while j <= limit {
sieve[j] = 1
j = j + pp
}
}
pp = pp + 1
}
free(sieve)
let mut i: i64 = 0
while i < pc {
let q: i64 = primes[i]
let mut offset: i64 = (q - (lo % q)) % q
let mut j: i64 = offset
while j < length {
is_comp[j] = 1
j = j + q
}
if lo <= q && q <= hi {
is_comp[q - lo] = 0
}
i = i + 1
}
free(primes)
let mut total: i64 = 0
let mut j: i64 = 0
while j < length {
if is_comp[j] == 0 {
let n: i64 = lo + j
if is_prime_64(n) != 0 {
total = (total + f_of_prime(n, modv)) % modv
}
}
j = j + 1
}
free(is_comp)
printf("%lld\n", total)
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 powmod_i64_i64_i64(int64_t a0, int64_t e0, int64_t m);
int32_t is_prime_64_i64(int64_t n);
int64_t gcd_u64_i64_i64(int64_t a0, int64_t b0);
int64_t rand64(void);
int64_t pollard_rho_i64(int64_t n);
int64_t factorize_i64_ptr_i64_ptr_i64(int64_t n0, int64_t* fp, int64_t* ep);
int64_t powmod_small_i64_i64_i64(int64_t a0, int64_t e0, int64_t m);
int64_t g_prime_power_mod_i64_i64_i64(int64_t q, int64_t e, int64_t modv);
int64_t g_from_factorization_ptr_i64_ptr_i64_i64_i64(int64_t* fp, int64_t* ep, int64_t nf, int64_t modv);
int64_t f_of_prime_i64_i64(int64_t p, int64_t modv);
int32_t main(void);
static const int64_t MOD = 993353399;
/* Module statics */
static int64_t rng_state = (-7046029254386353131);
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 powmod_i64_i64_i64(int64_t a0, int64_t e0, int64_t m) {
int64_t r = FLOW_CHECKED_MOD((1), (m));
int64_t a = FLOW_CHECKED_MOD((a0), (m));
int64_t e = e0;
while (e != 0) {
if (FLOW_CHECKED_MOD((e), (2)) == 1) {
r = mulmod_i64_i64_i64(r, a, m);
}
a = mulmod_i64_i64_i64(a, a, m);
e = FLOW_CHECKED_DIV((e), (2));
}
return r;
}
int32_t is_prime_64_i64(int64_t n) {
if (n < 2) {
return 0;
}
int64_t sp[12] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
int64_t i = 0;
while (i < 12) {
if (n == (((unsigned)(i) < 12) ? sp[i] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(i), 12), flow_fault_handler("array index out of bounds"), sp[0]))) {
return 1;
}
if (FLOW_CHECKED_MOD((n), ((((unsigned)(i) < 12) ? sp[i] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(i), 12), flow_fault_handler("array index out of bounds"), sp[0])))) == 0) {
return 0;
}
i = (i + 1);
}
int64_t d = (n - 1);
int64_t s = 0;
while (FLOW_CHECKED_MOD((d), (2)) == 0) {
d = FLOW_CHECKED_DIV((d), (2));
s = (s + 1);
}
int64_t bases[7] = { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 };
int64_t bi = 0;
while (bi < 7) {
int64_t a = (((unsigned)(bi) < 7) ? bases[bi] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(bi), 7), flow_fault_handler("array index out of bounds"), bases[0]));
if (FLOW_CHECKED_MOD((a), (n)) == 0) {
bi = (bi + 1);
} else {
int64_t x = powmod_i64_i64_i64(a, d, n);
if ((x == 1 || x == (n - 1))) {
bi = (bi + 1);
} else {
int32_t comp = 1;
int64_t j = 0;
while (j < (s - 1)) {
x = mulmod_i64_i64_i64(x, x, n);
if (x == (n - 1)) {
comp = 0;
j = s;
}
j = (j + 1);
}
if (comp != 0) {
return 0;
}
bi = (bi + 1);
}
}
}
return 1;
}
int64_t gcd_u64_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 rand64(void) {
int64_t x = rng_state;
x = (x ^ FLOW_CHECKED_SHR((x), (12)));
x = (x ^ FLOW_CHECKED_SHL((x), (25)));
x = (x ^ FLOW_CHECKED_SHR((x), (27)));
rng_state = x;
return (x * 2685821657736338717);
}
int64_t pollard_rho_i64(int64_t n) {
if (FLOW_CHECKED_MOD((n), (2)) == 0) {
return 2;
}
if (FLOW_CHECKED_MOD((n), (3)) == 0) {
return 3;
}
while (1) {
int64_t c = (FLOW_CHECKED_MOD((rand64()), ((n - 1))) + 1);
int64_t x = (FLOW_CHECKED_MOD((rand64()), ((n - 2))) + 2);
int64_t y = x;
int64_t d = 1;
while (d == 1) {
x = FLOW_CHECKED_MOD(((mulmod_i64_i64_i64(x, x, n) + c)), (n));
y = FLOW_CHECKED_MOD(((mulmod_i64_i64_i64(y, y, n) + c)), (n));
y = FLOW_CHECKED_MOD(((mulmod_i64_i64_i64(y, y, n) + c)), (n));
int64_t diff = (x - y);
if (diff < 0) {
diff = (-diff);
}
d = gcd_u64_i64_i64(diff, n);
}
if (d != n) {
return d;
}
}
return 0;
}
int64_t factorize_i64_ptr_i64_ptr_i64(int64_t n0, int64_t* fp, int64_t* ep) {
int64_t outn = 0;
int64_t n = n0;
if (n <= 1) {
return 0;
}
int64_t p = 2;
while (((p * p) <= n && p <= 10000)) {
if (FLOW_CHECKED_MOD((n), (p)) == 0) {
int64_t e = 0;
while (FLOW_CHECKED_MOD((n), (p)) == 0) {
n = FLOW_CHECKED_DIV((n), (p));
e = (e + 1);
}
fp[outn] = p;
ep[outn] = e;
outn = (outn + 1);
}
p = (p + 1);
}
if (n == 1) {
return outn;
}
if (is_prime_64_i64(n) != 0) {
fp[outn] = n;
ep[outn] = 1;
outn = (outn + 1);
return outn;
}
int64_t* stack = (int64_t*)(((int64_t*)(calloc(64, 8))));
int64_t sp = 0;
stack[sp] = n;
sp = (sp + 1);
while (sp > 0) {
sp = (sp - 1);
int64_t m = stack[sp];
if (m == 1) {
} else {
if (is_prime_64_i64(m) != 0) {
int32_t found = 0;
int64_t i = 0;
while (i < outn) {
if (fp[i] == m) {
ep[i] = (ep[i] + 1);
found = 1;
i = outn;
}
i = (i + 1);
}
if (found == 0) {
fp[outn] = m;
ep[outn] = 1;
outn = (outn + 1);
}
} else {
int64_t d = pollard_rho_i64(m);
stack[sp] = d;
sp = (sp + 1);
stack[sp] = FLOW_CHECKED_DIV((m), (d));
sp = (sp + 1);
}
}
}
free(stack);
return outn;
}
int64_t powmod_small_i64_i64_i64(int64_t a0, int64_t e0, int64_t m) {
int64_t r = FLOW_CHECKED_MOD((1), (m));
int64_t a = FLOW_CHECKED_MOD((a0), (m));
if (a < 0) {
a = (a + m);
}
int64_t e = e0;
while (e > 0) {
if (FLOW_CHECKED_MOD((e), (2)) == 1) {
r = FLOW_CHECKED_MOD(((r * a)), (m));
}
a = FLOW_CHECKED_MOD(((a * a)), (m));
e = FLOW_CHECKED_DIV((e), (2));
}
return r;
}
int64_t g_prime_power_mod_i64_i64_i64(int64_t q, int64_t e, int64_t modv) {
int64_t qm = FLOW_CHECKED_MOD((q), (modv));
int64_t s = 0;
int64_t t = 1;
while (t <= e) {
int64_t expv = (((3 * e) - t) - 2);
int64_t term = powmod_small_i64_i64_i64(qm, expv, modv);
s = FLOW_CHECKED_MOD(((s + (FLOW_CHECKED_MOD(((t * t)), (modv)) * term))), (modv));
t = (t + 1);
}
int64_t qm1 = FLOW_CHECKED_MOD(((q - 1)), (modv));
int64_t term1 = FLOW_CHECKED_MOD(((powmod_small_i64_i64_i64(qm1, 3, modv) * s)), (modv));
int64_t inner = FLOW_CHECKED_MOD((((e * (q - 1)) + q)), (modv));
int64_t term2 = FLOW_CHECKED_MOD(((powmod_small_i64_i64_i64(qm, ((2 * e) - 2), modv) * powmod_small_i64_i64_i64(inner, 2, modv))), (modv));
return FLOW_CHECKED_MOD(((term1 + term2)), (modv));
}
int64_t g_from_factorization_ptr_i64_ptr_i64_i64_i64(int64_t* fp, int64_t* ep, int64_t nf, int64_t modv) {
int64_t g = FLOW_CHECKED_MOD((1), (modv));
int64_t i = 0;
while (i < nf) {
g = FLOW_CHECKED_MOD(((g * g_prime_power_mod_i64_i64_i64(fp[i], ep[i], modv))), (modv));
i = (i + 1);
}
return g;
}
int64_t f_of_prime_i64_i64(int64_t p, int64_t modv) {
int64_t m = (p - 1);
int64_t* fp = (int64_t*)(((int64_t*)(calloc(64, 8))));
int64_t* ep = (int64_t*)(((int64_t*)(calloc(64, 8))));
int64_t nf = factorize_i64_ptr_i64_ptr_i64(m, fp, ep);
int64_t g = g_from_factorization_ptr_i64_ptr_i64_i64_i64(fp, ep, nf, modv);
free(fp);
free(ep);
int64_t mm = FLOW_CHECKED_MOD((m), (modv));
return FLOW_CHECKED_MOD((((mm * mm) + g)), (modv));
}
int32_t main(void) {
int64_t modv = MOD;
int64_t lo = 10000000000000000;
int64_t hi = (lo + 1000000);
int64_t length = ((hi - lo) + 1);
int8_t* is_comp = (int8_t*)(((int8_t*)(calloc(length, 1))));
int64_t limit = 200000;
int8_t* sieve = (int8_t*)(((int8_t*)(calloc((limit + 1), 1))));
int64_t* primes = (int64_t*)(((int64_t*)(calloc(20000, 8))));
int64_t pc = 0;
int64_t pp = 2;
while (pp <= limit) {
if (sieve[pp] == 0) {
primes[pc] = pp;
pc = (pc + 1);
int64_t j = (pp * pp);
while (j <= limit) {
sieve[j] = 1;
j = (j + pp);
}
}
pp = (pp + 1);
}
free(sieve);
int64_t i = 0;
while (i < pc) {
int64_t q = primes[i];
int64_t offset = FLOW_CHECKED_MOD(((q - FLOW_CHECKED_MOD((lo), (q)))), (q));
int64_t j = offset;
while (j < length) {
is_comp[j] = 1;
j = (j + q);
}
if ((lo <= q && q <= hi)) {
is_comp[(q - lo)] = 0;
}
i = (i + 1);
}
free(primes);
int64_t total = 0;
int64_t j = 0;
while (j < length) {
if (is_comp[j] == 0) {
int64_t n = (lo + j);
if (is_prime_64_i64(n) != 0) {
total = FLOW_CHECKED_MOD(((total + f_of_prime_i64_i64(n, modv))), (modv));
}
}
j = (j + 1);
}
free(is_comp);
printf("%lld\n", total);
return 0;
}