# Project Euler 238
# Infinite string tour: BBS digit stream, sum p(k) for k <= 2e15.
# Period digit-sum T ~ 8e7; cover residues with a u64 limb bitset.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
function pop64(x0: u64) -> i64 {
let mut x: u64 = x0
let mut c: i64 = 0
while x != 0 {
c = c + 1
x = x & (x - 1)
}
return c
}
function bit_set(bits: ptr<u64>, idx: i64) -> void {
let limb: i64 = idx >> 6
let bit: i64 = idx & 63
let one: u64 = 1
bits[limb] = bits[limb] | (one << (bit as u64))
}
function bits_copy(src: ptr<u64>, dst: ptr<u64>, nlimbs: i64) -> void {
let mut i: i64 = 0
while i < nlimbs {
dst[i] = src[i]
i = i + 1
}
}
function bits_zero(bits: ptr<u64>, nlimbs: i64) -> void {
let mut i: i64 = 0
let z: u64 = 0
while i < nlimbs {
bits[i] = z
i = i + 1
}
}
# Set every bit in [0, T).
function bits_fill(bits: ptr<u64>, T: i64) -> void {
let nlimbs: i64 = (T + 63) / 64
let all: u64 = 0
all = all - 1
let mut i: i64 = 0
let full: i64 = T / 64
while i < full {
bits[i] = all
i = i + 1
}
let remb: i64 = T & 63
if remb != 0 {
let one: u64 = 1
bits[full] = (one << (remb as u64)) - 1
}
}
# Circular right shift by d bits (1..9) of a T-bit set into dst.
function circ_rshift(src: ptr<u64>, dst: ptr<u64>, T: i64, d: i64) -> void {
let nlimbs: i64 = (T + 63) / 64
let du: u64 = d as u64
let one: u64 = 1
let low_mask: u64 = (one << du) - 1
let low: u64 = src[0] & low_mask
let sh: u64 = 64 - du
let mut i: i64 = 0
while i + 1 < nlimbs {
dst[i] = (src[i] >> du) | (src[i + 1] << sh)
i = i + 1
}
dst[nlimbs - 1] = src[nlimbs - 1] >> du
# Wrap the saved low d bits into positions [T-d, T).
let top: i64 = T - d
let limb: i64 = top >> 6
let bit: i64 = top & 63
let bu: u64 = bit as u64
dst[limb] = dst[limb] | (low << bu)
if bit + d > 64 {
dst[limb + 1] = dst[limb + 1] | (low >> (64 - bu))
}
# Mask to exactly T bits.
let remb: i64 = T & 63
if remb != 0 {
dst[nlimbs - 1] = dst[nlimbs - 1] & ((one << (remb as u64)) - 1)
}
}
# Popcount of bits of v0 whose global indices lie in [lo, hi], for limb li.
function pop_in_range(v0: u64, li: i64, lo: i64, hi: i64) -> i64 {
if hi < lo { return 0 }
let base: i64 = li * 64
if base > hi { return 0 }
if base + 63 < lo { return 0 }
let mut start: i64 = lo - base
if start < 0 { start = 0 }
let mut endb: i64 = hi - base
if endb > 63 { endb = 63 }
let one: u64 = 1
let width: i64 = endb - start + 1
let mut mask: u64 = 0
if width == 64 {
mask = mask - 1
} else {
mask = ((one << (width as u64)) - 1) << (start as u64)
}
return pop64(v0 & mask)
}
# Write decimal digits of s0 at digits[pos0..]; return new pos. Digit sum via sum_out.
function digit_sum_and_emit(s0: i64, digits: ptr<i8>, pos0: i64, sum_out: ptr<i64>) -> i64 {
let mut tmpd: array<i32, 8> = [0, 0, 0, 0, 0, 0, 0, 0]
let mut t: i64 = s0
let mut len: i32 = 0
if t == 0 {
digits[pos0] = 0
return pos0 + 1
}
while t > 0 {
tmpd[len] = (t % 10) as i32
t = t / 10
len = len + 1
}
let mut pos: i64 = pos0
let mut di: i32 = len - 1
while di >= 0 {
let d: i32 = tmpd[di]
digits[pos] = d as i8
sum_out[0] = sum_out[0] + (d as i64)
pos = pos + 1
di = di - 1
}
return pos
}
function bbs_period(s0: i64, mod: i64) -> i64 {
let mut s: i64 = s0
let mut n: i64 = 0
while n < 100000000 {
s = (s * s) % mod
n = n + 1
if s == s0 { return n }
}
return 0
}
function main() -> i32 {
let S0: i64 = 14025256
let MOD: i64 = 20300713
let TARGET: i64 = 2000000000000000
let period: i64 = bbs_period(S0, MOD)
# One period has < 8 digits per term.
let dig_cap: i64 = period * 8 + 8
let digits: ptr<i8> = calloc(dig_cap, 1)
if digits == null { return 1 }
let mut s: i64 = S0
let mut L: i64 = 0
let mut Tacc: i64 = 0
let mut i: i64 = 0
while i < period {
L = digit_sum_and_emit(s, digits, L, &Tacc)
s = (s * s) % MOD
i = i + 1
}
let T: i64 = Tacc
let nlimbs: i64 = (T + 63) / 64
let present: ptr<u64> = calloc(nlimbs, 8)
let unknown: ptr<u64> = calloc(nlimbs, 8)
let rot: ptr<u64> = calloc(nlimbs, 8)
let tmp: ptr<u64> = calloc(nlimbs, 8)
if present == null || unknown == null || rot == null || tmp == null { return 1 }
bits_zero(present, nlimbs)
bit_set(present, 0)
let mut ps: i64 = 0
i = 0
while i < L {
ps = ps + (digits[i] as i64)
if ps == T {
bit_set(present, 0)
} else {
bit_set(present, ps)
}
i = i + 1
}
bits_fill(unknown, T)
bits_copy(present, rot, nlimbs)
let rem: i64 = TARGET % T
let mut sum_all: i64 = 0
let mut sum_1000: i64 = 0
let mut sum_rem: i64 = 0
let mut left: i64 = T
let mut z: i64 = 0
while z < L {
let mut cnt: i64 = 0
let mut c1000: i64 = 0
let mut crem: i64 = 0
let mut li: i64 = 0
while li < nlimbs {
let v: u64 = unknown[li] & rot[li]
if v != 0 {
cnt = cnt + pop64(v)
c1000 = c1000 + pop_in_range(v, li, 1, 1000)
crem = crem + pop_in_range(v, li, 1, rem)
unknown[li] = unknown[li] ^ v
}
li = li + 1
}
let w: i64 = z + 1
if cnt != 0 {
sum_all = sum_all + w * cnt
sum_1000 = sum_1000 + w * c1000
sum_rem = sum_rem + w * crem
left = left - cnt
if left == 0 { break }
}
let d: i64 = digits[z] as i64
if d != 0 {
bits_zero(tmp, nlimbs)
circ_rshift(rot, tmp, T, d)
bits_copy(tmp, rot, nlimbs)
}
z = z + 1
}
if sum_1000 != 4742 {
printf("check failed: %lld\n", sum_1000)
return 1
}
let q: i64 = TARGET / T
let ans: i64 = q * sum_all + sum_rem
printf("%lld\n", ans)
free(digits)
free(present)
free(unknown)
free(rot)
free(tmp)
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 pop64_u64(uint64_t x0);
void bit_set_ptr_u64_i64(uint64_t* bits, int64_t idx);
void bits_copy_ptr_u64_ptr_u64_i64(uint64_t* src, uint64_t* dst, int64_t nlimbs);
void bits_zero_ptr_u64_i64(uint64_t* bits, int64_t nlimbs);
void bits_fill_ptr_u64_i64(uint64_t* bits, int64_t T);
void circ_rshift_ptr_u64_ptr_u64_i64_i64(uint64_t* src, uint64_t* dst, int64_t T, int64_t d);
int64_t pop_in_range_u64_i64_i64_i64(uint64_t v0, int64_t li, int64_t lo, int64_t hi);
int64_t digit_sum_and_emit_i64_ptr_i8_i64_ptr_i64(int64_t s0, int8_t* digits, int64_t pos0, int64_t* sum_out);
int64_t bbs_period_i64_i64(int64_t s0, int64_t mod);
int32_t main(void);
int64_t pop64_u64(uint64_t x0) {
uint64_t x = x0;
int64_t c = 0;
while (x != 0) {
c = (c + 1);
x = (x & (x - 1));
}
return c;
}
void bit_set_ptr_u64_i64(uint64_t* bits, int64_t idx) {
int64_t limb = FLOW_CHECKED_SHR((idx), (6));
int64_t bit = (idx & 63);
uint64_t one = 1;
bits[limb] = (bits[limb] | FLOW_CHECKED_SHL((one), (((uint64_t)(bit)))));
}
void bits_copy_ptr_u64_ptr_u64_i64(uint64_t* src, uint64_t* dst, int64_t nlimbs) {
int64_t i = 0;
while (i < nlimbs) {
dst[i] = src[i];
i = (i + 1);
}
}
void bits_zero_ptr_u64_i64(uint64_t* bits, int64_t nlimbs) {
int64_t i = 0;
uint64_t z = 0;
while (i < nlimbs) {
bits[i] = z;
i = (i + 1);
}
}
void bits_fill_ptr_u64_i64(uint64_t* bits, int64_t T) {
int64_t nlimbs = FLOW_CHECKED_DIV(((T + 63)), (64));
uint64_t all = 0;
all = (all - 1);
int64_t i = 0;
int64_t full = FLOW_CHECKED_DIV((T), (64));
while (i < full) {
bits[i] = all;
i = (i + 1);
}
int64_t remb = (T & 63);
if (remb != 0) {
uint64_t one = 1;
bits[full] = (FLOW_CHECKED_SHL((one), (((uint64_t)(remb)))) - 1);
}
}
void circ_rshift_ptr_u64_ptr_u64_i64_i64(uint64_t* src, uint64_t* dst, int64_t T, int64_t d) {
int64_t nlimbs = FLOW_CHECKED_DIV(((T + 63)), (64));
uint64_t du = ((uint64_t)(d));
uint64_t one = 1;
uint64_t low_mask = (FLOW_CHECKED_SHL((one), (du)) - 1);
uint64_t low = (src[0] & low_mask);
uint64_t sh = (64 - du);
int64_t i = 0;
while ((i + 1) < nlimbs) {
dst[i] = (FLOW_CHECKED_SHR((src[i]), (du)) | FLOW_CHECKED_SHL((src[(i + 1)]), (sh)));
i = (i + 1);
}
dst[(nlimbs - 1)] = FLOW_CHECKED_SHR((src[(nlimbs - 1)]), (du));
int64_t top = (T - d);
int64_t limb = FLOW_CHECKED_SHR((top), (6));
int64_t bit = (top & 63);
uint64_t bu = ((uint64_t)(bit));
dst[limb] = (dst[limb] | FLOW_CHECKED_SHL((low), (bu)));
if ((bit + d) > 64) {
dst[(limb + 1)] = (dst[(limb + 1)] | FLOW_CHECKED_SHR((low), ((64 - bu))));
}
int64_t remb = (T & 63);
if (remb != 0) {
dst[(nlimbs - 1)] = (dst[(nlimbs - 1)] & (FLOW_CHECKED_SHL((one), (((uint64_t)(remb)))) - 1));
}
}
int64_t pop_in_range_u64_i64_i64_i64(uint64_t v0, int64_t li, int64_t lo, int64_t hi) {
if (hi < lo) {
return 0;
}
int64_t base = (li * 64);
if (base > hi) {
return 0;
}
if ((base + 63) < lo) {
return 0;
}
int64_t start = (lo - base);
if (start < 0) {
start = 0;
}
int64_t endb = (hi - base);
if (endb > 63) {
endb = 63;
}
uint64_t one = 1;
int64_t width = ((endb - start) + 1);
uint64_t mask = 0;
if (width == 64) {
mask = (mask - 1);
} else {
mask = FLOW_CHECKED_SHL(((FLOW_CHECKED_SHL((one), (((uint64_t)(width)))) - 1)), (((uint64_t)(start))));
}
return pop64_u64((v0 & mask));
}
int64_t digit_sum_and_emit_i64_ptr_i8_i64_ptr_i64(int64_t s0, int8_t* digits, int64_t pos0, int64_t* sum_out) {
int32_t tmpd[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
int64_t t = s0;
int32_t len = 0;
if (t == 0) {
digits[pos0] = 0;
return (pos0 + 1);
}
while (t > 0) {
tmpd[len] = ((int32_t)(FLOW_CHECKED_MOD((t), (10))));
t = FLOW_CHECKED_DIV((t), (10));
len = (len + 1);
}
int64_t pos = pos0;
int32_t di = (len - 1);
while (di >= 0) {
int32_t d = (((unsigned)(di) < 8) ? tmpd[di] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(di), 8), flow_fault_handler("array index out of bounds"), tmpd[0]));
digits[pos] = ((int8_t)(d));
sum_out[0] = (sum_out[0] + ((int64_t)(d)));
pos = (pos + 1);
di = (di - 1);
}
return pos;
}
int64_t bbs_period_i64_i64(int64_t s0, int64_t mod) {
int64_t s = s0;
int64_t n = 0;
while (n < 100000000) {
s = FLOW_CHECKED_MOD(((s * s)), (mod));
n = (n + 1);
if (s == s0) {
return n;
}
}
return 0;
}
int32_t main(void) {
int64_t S0 = 14025256;
int64_t MOD = 20300713;
int64_t TARGET = 2000000000000000;
int64_t period = bbs_period_i64_i64(S0, MOD);
int64_t dig_cap = ((period * 8) + 8);
int8_t* digits = (int8_t*)(calloc(dig_cap, 1));
if (digits == NULL) {
return 1;
}
int64_t s = S0;
int64_t L = 0;
int64_t Tacc = 0;
int64_t i = 0;
while (i < period) {
L = digit_sum_and_emit_i64_ptr_i8_i64_ptr_i64(s, digits, L, (&(Tacc)));
s = FLOW_CHECKED_MOD(((s * s)), (MOD));
i = (i + 1);
}
int64_t T = Tacc;
int64_t nlimbs = FLOW_CHECKED_DIV(((T + 63)), (64));
uint64_t* present = (uint64_t*)(calloc(nlimbs, 8));
uint64_t* unknown = (uint64_t*)(calloc(nlimbs, 8));
uint64_t* rot = (uint64_t*)(calloc(nlimbs, 8));
uint64_t* tmp = (uint64_t*)(calloc(nlimbs, 8));
if ((((present == NULL || unknown == NULL) || rot == NULL) || tmp == NULL)) {
return 1;
}
bits_zero_ptr_u64_i64(present, nlimbs);
bit_set_ptr_u64_i64(present, 0);
int64_t ps = 0;
i = 0;
while (i < L) {
ps = (ps + ((int64_t)(digits[i])));
if (ps == T) {
bit_set_ptr_u64_i64(present, 0);
} else {
bit_set_ptr_u64_i64(present, ps);
}
i = (i + 1);
}
bits_fill_ptr_u64_i64(unknown, T);
bits_copy_ptr_u64_ptr_u64_i64(present, rot, nlimbs);
int64_t rem = FLOW_CHECKED_MOD((TARGET), (T));
int64_t sum_all = 0;
int64_t sum_1000 = 0;
int64_t sum_rem = 0;
int64_t left = T;
int64_t z = 0;
while (z < L) {
int64_t cnt = 0;
int64_t c1000 = 0;
int64_t crem = 0;
int64_t li = 0;
while (li < nlimbs) {
uint64_t v = (unknown[li] & rot[li]);
if (v != 0) {
cnt = (cnt + pop64_u64(v));
c1000 = (c1000 + pop_in_range_u64_i64_i64_i64(v, li, 1, 1000));
crem = (crem + pop_in_range_u64_i64_i64_i64(v, li, 1, rem));
unknown[li] = (unknown[li] ^ v);
}
li = (li + 1);
}
int64_t w = (z + 1);
if (cnt != 0) {
sum_all = (sum_all + (w * cnt));
sum_1000 = (sum_1000 + (w * c1000));
sum_rem = (sum_rem + (w * crem));
left = (left - cnt);
if (left == 0) {
break;
}
}
int64_t d = ((int64_t)(digits[z]));
if (d != 0) {
bits_zero_ptr_u64_i64(tmp, nlimbs);
circ_rshift_ptr_u64_ptr_u64_i64_i64(rot, tmp, T, d);
bits_copy_ptr_u64_ptr_u64_i64(tmp, rot, nlimbs);
}
z = (z + 1);
}
if (sum_1000 != 4742) {
printf("check failed: %lld\n", sum_1000);
return 1;
}
int64_t q = FLOW_CHECKED_DIV((TARGET), (T));
int64_t ans = ((q * sum_all) + sum_rem);
printf("%lld\n", ans);
free(digits);
free(present);
free(unknown);
free(rot);
free(tmp);
return 0;
}