# Project Euler 969
# Sum S(n) for n=1..10^18 via Lagrange interpolation.
# Ported from native C to pure Flow.
import euler.nt { mod_pow }
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>)
function printf(fmt: ptr<i8>, ...) -> i32
}
const MOD: i64 = 1000000007
let mut primes: ptr<i32> = null
let mut num_primes: i32 = 0
function sieve(n: i32) -> void {
num_primes = 0
let is_composite: ptr<i8> = calloc((n + 1) as i64, 1) as ptr<i8>
let mut i: i32 = 2
while i <= n {
if is_composite[i] == 0 {
num_primes = num_primes + 1
let mut j: i32 = i * i
while j <= n {
is_composite[j] = 1
j = j + i
}
}
i = i + 1
}
primes = calloc(num_primes as i64, 4) as ptr<i32>
let mut k: i32 = 0
i = 2
while i <= n {
if is_composite[i] == 0 {
primes[k] = i
k = k + 1
}
i = i + 1
}
free(is_composite as ptr<void>)
}
function modinv(x: i64) -> i64 {
return mod_pow(x % MOD, MOD - 2, MOD)
}
# sum_{t=1}^T t^k mod MOD using Lagrange interpolation
function sum_pows_lagrange(T: i64, k: i32) -> i64 {
let d: i32 = k + 1
if T <= (d as i64) {
let mut s: i64 = 0
let mut i: i64 = 1
while i <= T {
s = (s + mod_pow(i, k as i64, MOD)) % MOD
i = i + 1
}
return s
}
# y_i = sum_{t=1}^i t^k for i=0..d
let ys: ptr<i64> = calloc(60, 8) as ptr<i64>
ys[0] = 0
let mut i2: i32 = 1
while i2 <= d {
ys[i2] = (ys[i2 - 1] + mod_pow(i2 as i64, k as i64, MOD)) % MOD
i2 = i2 + 1
}
let fact: ptr<i64> = calloc(60, 8) as ptr<i64>
fact[0] = 1
let mut i3: i32 = 1
while i3 <= d {
fact[i3] = fact[i3 - 1] * (i3 as i64) % MOD
i3 = i3 + 1
}
let invfact: ptr<i64> = calloc(60, 8) as ptr<i64>
invfact[d] = modinv(fact[d])
let mut i4: i32 = d
while i4 >= 1 {
invfact[i4 - 1] = invfact[i4] * (i4 as i64) % MOD
i4 = i4 - 1
}
let Tmod: i64 = T % MOD
# pre[i] = prod_{j=0..i-1} (T - j)
let pre: ptr<i64> = calloc(60, 8) as ptr<i64>
pre[0] = 1
let mut i5: i32 = 1
while i5 <= d {
let val: i64 = (Tmod - ((i5 - 1) as i64)) % MOD
pre[i5] = pre[i5 - 1] * val % MOD
if pre[i5] < 0 { pre[i5] = pre[i5] + MOD }
i5 = i5 + 1
}
# suf[i] = prod_{j=i..d} (T - j)
let suf: ptr<i64> = calloc(62, 8) as ptr<i64>
suf[d + 1] = 1
let mut i6: i32 = d
while i6 >= 0 {
let val: i64 = (Tmod - (i6 as i64)) % MOD
suf[i6] = suf[i6 + 1] * val % MOD
if suf[i6] < 0 { suf[i6] = suf[i6] + MOD }
i6 = i6 - 1
}
let mut res: i64 = 0
let mut i7: i32 = 0
while i7 <= d {
let num: i64 = pre[i7] * suf[i7 + 1] % MOD
let mut denom: i64 = fact[i7] * fact[d - i7] % MOD
if ((d - i7) & 1) != 0 {
denom = (MOD - denom) % MOD
}
let invden: i64 = modinv(denom)
let li: i64 = num * invden % MOD
res = (res + ys[i7] * li) % MOD
i7 = i7 + 1
}
free(ys as ptr<void>)
free(fact as ptr<void>)
free(invfact as ptr<void>)
free(pre as ptr<void>)
free(suf as ptr<void>)
return res
}
# Compute M_k = prod_{p <= k} p^{ceil(v_p(k!)/k)}
function compute_Mk(k: i32) -> i128 {
if k == 0 { return 1 }
let mut M: i128 = 1
let mut pi: i32 = 0
while pi < num_primes {
let p: i32 = primes[pi]
if p > k { break }
let mut v: i32 = 0
let mut pp: i128 = p as i128
while pp <= (k as i128) {
v = v + (k as i32) / (pp as i32)
pp = pp * (p as i128)
}
let e: i32 = (v + k - 1) / k
if e > 0 {
let mut pe: i128 = 1
let mut i: i32 = 0
while i < e {
pe = pe * (p as i128)
i = i + 1
}
M = M * pe
}
pi = pi + 1
}
return M
}
function compute_sum_S_upto_N(N: i64) -> i64 {
let mut total: i64 = 0
# Precompute factorials mod MOD up to 200
let MAXF: i32 = 200
let fact: ptr<i64> = calloc(201, 8) as ptr<i64>
fact[0] = 1
let mut i: i32 = 1
while i <= MAXF {
fact[i] = fact[i - 1] * (i as i64) % MOD
i = i + 1
}
let mut k: i32 = 0
while true {
let M: i128 = compute_Mk(k)
if M > (N as i128) { break }
let T: i64
if k == 0 {
T = N
} else {
if N < (k as i64) { break }
T = (N - (k as i64)) / (M as i64)
if T <= 0 {
k = k + 1
continue
}
}
# C_k = (-1)^k * M^k / k! mod MOD
let C: i64
if k == 0 {
C = 1
} else {
C = mod_pow((M % (MOD as i128)) as i64, k as i64, MOD) * modinv(fact[k]) % MOD
if (k & 1) != 0 {
C = (MOD - C) % MOD
}
}
let sum_tk: i64 = sum_pows_lagrange(T, k)
let contrib: i64 = C * sum_tk % MOD
total = (total + contrib) % MOD
k = k + 1
}
free(fact as ptr<void>)
return total
}
function main() -> i32 {
sieve(1000)
let result: i64 = compute_sum_S_upto_N(1000000000000000000)
printf("%lld\n", result)
free(primes as ptr<void>)
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 lcm_i64_i64(int64_t a, int64_t b);
int64_t isqrt_i64(int64_t n);
int64_t mulmod_i64_i64_i64(int64_t a0, int64_t b0, int64_t mod);
int64_t mod_pow_i64_i64_i64(int64_t base, int64_t exp, int64_t mod);
bool is_prime_i64(int64_t n);
void sieve_i32(int32_t n);
int64_t modinv_i64(int64_t x);
int64_t sum_pows_lagrange_i64_i32(int64_t T, int32_t k);
__int128 compute_Mk_i32(int32_t k);
int64_t compute_sum_S_upto_N_i64(int64_t N);
int32_t main(void);
static const int64_t MOD = 1000000007;
/* Module statics */
static int32_t* primes = NULL;
static int32_t num_primes = 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 lcm_i64_i64(int64_t a, int64_t b) {
if ((a == 0 || b == 0)) {
return 0;
}
return (FLOW_CHECKED_DIV((a), (gcd_i64_i64(a, b))) * b);
}
int64_t isqrt_i64(int64_t n) {
if (n < 2) {
return n;
}
int64_t x = n;
int64_t y = FLOW_CHECKED_DIV(((x + 1)), (2));
while (y < x) {
x = y;
y = FLOW_CHECKED_DIV(((x + FLOW_CHECKED_DIV((n), (x)))), (2));
}
return x;
}
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));
int64_t result = 0;
while (b > 0) {
if (FLOW_CHECKED_MOD((b), (2)) == 1) {
result = FLOW_CHECKED_MOD(((result + a)), (mod));
}
a = FLOW_CHECKED_MOD(((a * 2)), (mod));
b = FLOW_CHECKED_DIV((b), (2));
}
return result;
}
int64_t mod_pow_i64_i64_i64(int64_t base, int64_t exp, int64_t mod) {
if (mod == 1) {
return 0;
}
int64_t result = 1;
int64_t b = FLOW_CHECKED_MOD((base), (mod));
int64_t e = exp;
while (e > 0) {
if (FLOW_CHECKED_MOD((e), (2)) == 1) {
result = mulmod_i64_i64_i64(result, b, mod);
}
b = mulmod_i64_i64_i64(b, b, mod);
e = FLOW_CHECKED_DIV((e), (2));
}
return result;
}
bool is_prime_i64(int64_t n) {
if (n < 2) {
return 0;
}
if (n < 4) {
return 1;
}
if ((FLOW_CHECKED_MOD((n), (2)) == 0 || FLOW_CHECKED_MOD((n), (3)) == 0)) {
return 0;
}
int64_t i = 5;
while ((i * i) <= n) {
if ((FLOW_CHECKED_MOD((n), (i)) == 0 || FLOW_CHECKED_MOD((n), ((i + 2))) == 0)) {
return 0;
}
i = (i + 6);
}
return 1;
}
void sieve_i32(int32_t n) {
num_primes = 0;
int8_t* is_composite = (int8_t*)(((int8_t*)(calloc(((int64_t)((n + 1))), 1))));
int32_t i = 2;
while (i <= n) {
if (is_composite[i] == 0) {
num_primes = (num_primes + 1);
int32_t j = (i * i);
while (j <= n) {
is_composite[j] = 1;
j = (j + i);
}
}
i = (i + 1);
}
primes = ((int32_t*)(calloc(((int64_t)(num_primes)), 4)));
int32_t k = 0;
i = 2;
while (i <= n) {
if (is_composite[i] == 0) {
primes[k] = i;
k = (k + 1);
}
i = (i + 1);
}
free(((void*)(is_composite)));
}
int64_t modinv_i64(int64_t x) {
return mod_pow_i64_i64_i64(FLOW_CHECKED_MOD((x), (MOD)), (MOD - 2), MOD);
}
int64_t sum_pows_lagrange_i64_i32(int64_t T, int32_t k) {
int32_t d = (k + 1);
if (T <= ((int64_t)(d))) {
int64_t s = 0;
int64_t i = 1;
while (i <= T) {
s = FLOW_CHECKED_MOD(((s + mod_pow_i64_i64_i64(i, ((int64_t)(k)), MOD))), (MOD));
i = (i + 1);
}
return s;
}
int64_t* ys = (int64_t*)(((int64_t*)(calloc(60, 8))));
ys[0] = 0;
int32_t i2 = 1;
while (i2 <= d) {
ys[i2] = FLOW_CHECKED_MOD(((ys[(i2 - 1)] + mod_pow_i64_i64_i64(((int64_t)(i2)), ((int64_t)(k)), MOD))), (MOD));
i2 = (i2 + 1);
}
int64_t* fact = (int64_t*)(((int64_t*)(calloc(60, 8))));
fact[0] = 1;
int32_t i3 = 1;
while (i3 <= d) {
fact[i3] = FLOW_CHECKED_MOD(((fact[(i3 - 1)] * ((int64_t)(i3)))), (MOD));
i3 = (i3 + 1);
}
int64_t* invfact = (int64_t*)(((int64_t*)(calloc(60, 8))));
invfact[d] = modinv_i64(fact[d]);
int32_t i4 = d;
while (i4 >= 1) {
invfact[(i4 - 1)] = FLOW_CHECKED_MOD(((invfact[i4] * ((int64_t)(i4)))), (MOD));
i4 = (i4 - 1);
}
int64_t Tmod = FLOW_CHECKED_MOD((T), (MOD));
int64_t* pre = (int64_t*)(((int64_t*)(calloc(60, 8))));
pre[0] = 1;
int32_t i5 = 1;
while (i5 <= d) {
int64_t val = FLOW_CHECKED_MOD(((Tmod - ((int64_t)((i5 - 1))))), (MOD));
pre[i5] = FLOW_CHECKED_MOD(((pre[(i5 - 1)] * val)), (MOD));
if (pre[i5] < 0) {
pre[i5] = (pre[i5] + MOD);
}
i5 = (i5 + 1);
}
int64_t* suf = (int64_t*)(((int64_t*)(calloc(62, 8))));
suf[(d + 1)] = 1;
int32_t i6 = d;
while (i6 >= 0) {
int64_t val = FLOW_CHECKED_MOD(((Tmod - ((int64_t)(i6)))), (MOD));
suf[i6] = FLOW_CHECKED_MOD(((suf[(i6 + 1)] * val)), (MOD));
if (suf[i6] < 0) {
suf[i6] = (suf[i6] + MOD);
}
i6 = (i6 - 1);
}
int64_t res = 0;
int32_t i7 = 0;
while (i7 <= d) {
int64_t num = FLOW_CHECKED_MOD(((pre[i7] * suf[(i7 + 1)])), (MOD));
int64_t denom = FLOW_CHECKED_MOD(((fact[i7] * fact[(d - i7)])), (MOD));
if (((d - i7) & 1) != 0) {
denom = FLOW_CHECKED_MOD(((MOD - denom)), (MOD));
}
int64_t invden = modinv_i64(denom);
int64_t li = FLOW_CHECKED_MOD(((num * invden)), (MOD));
res = FLOW_CHECKED_MOD(((res + (ys[i7] * li))), (MOD));
i7 = (i7 + 1);
}
free(((void*)(ys)));
free(((void*)(fact)));
free(((void*)(invfact)));
free(((void*)(pre)));
free(((void*)(suf)));
return res;
}
__int128 compute_Mk_i32(int32_t k) {
if (k == 0) {
return 1;
}
__int128 M = 1;
int32_t pi = 0;
while (pi < num_primes) {
int32_t p = primes[pi];
if (p > k) {
break;
}
int32_t v = 0;
__int128 pp = ((__int128)(p));
while (pp <= ((__int128)(k))) {
v = (v + FLOW_CHECKED_DIV((((int32_t)(k))), (((int32_t)(pp)))));
pp = (pp * ((__int128)(p)));
}
int32_t e = FLOW_CHECKED_DIV((((v + k) - 1)), (k));
if (e > 0) {
__int128 pe = 1;
int32_t i = 0;
while (i < e) {
pe = (pe * ((__int128)(p)));
i = (i + 1);
}
M = (M * pe);
}
pi = (pi + 1);
}
return M;
}
int64_t compute_sum_S_upto_N_i64(int64_t N) {
int64_t total = 0;
int32_t MAXF = 200;
int64_t* fact = (int64_t*)(((int64_t*)(calloc(201, 8))));
fact[0] = 1;
int32_t i = 1;
while (i <= MAXF) {
fact[i] = FLOW_CHECKED_MOD(((fact[(i - 1)] * ((int64_t)(i)))), (MOD));
i = (i + 1);
}
int32_t k = 0;
while (1) {
__int128 M = compute_Mk_i32(k);
if (M > ((__int128)(N))) {
break;
}
int64_t T;
if (k == 0) {
T = N;
} else {
if (N < ((int64_t)(k))) {
break;
}
T = FLOW_CHECKED_DIV(((N - ((int64_t)(k)))), (((int64_t)(M))));
if (T <= 0) {
k = (k + 1);
continue;
}
}
int64_t C;
if (k == 0) {
C = 1;
} else {
C = FLOW_CHECKED_MOD(((mod_pow_i64_i64_i64(((int64_t)(FLOW_CHECKED_MOD((M), (((__int128)(MOD)))))), ((int64_t)(k)), MOD) * modinv_i64(fact[k]))), (MOD));
if ((k & 1) != 0) {
C = FLOW_CHECKED_MOD(((MOD - C)), (MOD));
}
}
int64_t sum_tk = sum_pows_lagrange_i64_i32(T, k);
int64_t contrib = FLOW_CHECKED_MOD(((C * sum_tk)), (MOD));
total = FLOW_CHECKED_MOD(((total + contrib)), (MOD));
k = (k + 1);
}
free(((void*)(fact)));
return total;
}
int32_t main(void) {
sieve_i32(1000);
int64_t result = compute_sum_S_upto_N_i64(1000000000000000000);
printf("%lld\n", result);
free(((void*)(primes)));
return 0;
}