Sum of B(n^2) for n=1..10^16-1, mod 1e9+7. B(n) = next lexicographic permutation of n's digits (0 if none). Pure Flow port of the native C solver using i128.
# Project Euler 925: Next permutation of squares
# Sum of B(n^2) for n=1..10^16-1, mod 1e9+7.
# B(n) = next lexicographic permutation of n's digits (0 if none).
# Pure Flow port of the native C solver using i128.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
function malloc(n: i64) -> ptr<void>
function printf(fmt: ptr<i8>, ...) -> i32
}
const MOD: i64 = 1000000007
let mut g_pow10: ptr<i128> = null as ptr<i128>
let mut g_pow10_mod: ptr<i64> = null as ptr<i64>
let mut g_K: i32 = 0
# Compute (next_permutation(n) - n) % MOD for n given as i128.
function delta_mod_i128(n0: i128) -> i64 {
if n0 == 0 { return 0 }
let digits: array<i32, 40> = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
let mut len: i32 = 0
let mut tmp: i128 = n0
while tmp > 0 {
digits[len] = (tmp % 10) as i32
len = len + 1
tmp = tmp / 10
}
# Reverse to most-significant first
let mut i: i32 = 0
let mut j: i32 = len - 1
while i < j {
let t: i32 = digits[i]
digits[i] = digits[j]
digits[j] = t
i = i + 1
j = j - 1
}
# n mod MOD from digits
let mut n_mod: i64 = 0
for d in 0..len {
n_mod = ((n_mod * 10) % MOD + (digits[d] as i64)) % MOD
}
# Find largest i with digits[i] < digits[i+1]
let mut ii: i32 = len - 2
while ii >= 0 && digits[ii] >= digits[ii + 1] {
ii = ii - 1
}
if ii < 0 {
let r: i64 = (MOD - n_mod) % MOD
return r
}
# Find largest j > i with digits[j] > digits[i]
let mut jj: i32 = len - 1
while digits[jj] <= digits[ii] {
jj = jj - 1
}
# Swap
let t: i32 = digits[ii]
digits[ii] = digits[jj]
digits[jj] = t
# Reverse suffix starting at i+1
let mut lo: i32 = ii + 1
let mut hi: i32 = len - 1
while lo < hi {
let tt: i32 = digits[lo]
digits[lo] = digits[hi]
digits[hi] = tt
lo = lo + 1
hi = hi - 1
}
# next_permutation value mod MOD
let mut np_mod: i64 = 0
for d in 0..len {
np_mod = ((np_mod * 10) % MOD + (digits[d] as i64)) % MOD
}
let dm: i64 = (np_mod - n_mod) % MOD
if dm < 0 { return dm + MOD }
return dm
}
function mod_pow(a0: i64, e0: i64) -> i64 {
let mut r: i64 = 1
let mut a: i64 = a0 % MOD
if a < 0 { a = a + MOD }
let mut e: i64 = e0
while e > 0 {
if (e & 1) != 0 {
r = ((r * a) % MOD) as i64
}
a = ((a * a) % MOD) as i64
e = e >> 1
}
return r
}
function square_sum_below_power10(k: i32) -> i64 {
let n: i64 = mod_pow(10, (k as i64))
let mut t1: i64 = (n - 1) % MOD
if t1 < 0 { t1 = t1 + MOD }
let mut t2: i64 = (2 * n - 1) % MOD
if t2 < 0 { t2 = t2 + MOD }
let inv6: i64 = mod_pow(6, MOD - 2)
let mut result: i64 = ((n * t1) % MOD) as i64
result = ((result * t2) % MOD) as i64
result = ((result * inv6) % MOD) as i64
return result
}
function recurse(suffix0: i64, width: i32, trailing_zeros: i32) -> i64 {
let k: i32 = g_K
if width + trailing_zeros >= k {
let val: i128 = (suffix0 as i128) * g_pow10[trailing_zeros]
let sq: i128 = val * val
return delta_mod_i128(sq)
}
let square_suffix: i64 = ((suffix0 as i128) * (suffix0 as i128) % g_pow10[width]) as i64
let left_digit: i64 = square_suffix / (g_pow10[width - 1] as i64)
let next_place: i128 = g_pow10[width]
let next_modulus: i128 = g_pow10[width + 1]
let free_digits: i32 = k - width - trailing_zeros - 1
let completion_count: i64 = g_pow10_mod[free_digits]
let mut total: i64 = 0
for digit in 0..10 {
let next_suffix: i128 = (digit as i128) * next_place + (suffix0 as i128)
let next_sq: i128 = next_suffix * next_suffix
let new_digit: i64 = ((next_sq % next_modulus) / next_place) as i64
if new_digit >= left_digit {
total = (total + recurse((next_suffix as i64), width + 1, trailing_zeros)) % MOD
} else {
let next_sq_mod: i128 = next_sq % next_modulus
let visible: i128 = next_sq_mod * g_pow10[2 * trailing_zeros]
let representative: i128 = g_pow10[width + 1 + 2 * trailing_zeros] + visible
let representative_delta: i64 = delta_mod_i128(representative)
if next_sq < next_modulus {
total = (total + delta_mod_i128(visible)) % MOD
let mut term: i64 = (completion_count - 1) % MOD
if term < 0 { term = term + MOD }
term = ((term * representative_delta) % MOD) as i64
total = (total + term) % MOD
} else {
let term: i64 = ((completion_count * representative_delta) % MOD) as i64
total = (total + term) % MOD
}
}
}
return total
}
function correction_sum(k: i32) -> i64 {
g_K = k
# Build pow10 array (exact, up to 2*k+2) using i128
g_pow10[0] = 1
for i in 1..(2 * k + 3) {
g_pow10[i] = g_pow10[i - 1] * 10
}
# Build pow10_mod (up to k)
g_pow10_mod[0] = 1
for i in 1..(k + 1) {
g_pow10_mod[i] = (g_pow10_mod[i - 1] * 10) % MOD
}
let mut total: i64 = 0
for trailing_zeros in 0..k {
for last_digit in 1..10 {
total = (total + recurse((last_digit as i64), 1, trailing_zeros)) % MOD
}
}
return total
}
function main() -> i32 {
g_pow10 = (malloc(40 * 16)) as ptr<i128>
g_pow10_mod = (malloc(20 * 8)) as ptr<i64>
let k: i32 = 16
let result: i64 = (square_sum_below_power10(k) + correction_sum(k)) % MOD
printf("%lld\n", result)
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 delta_mod_i128_i128(__int128 n0);
int64_t mod_pow_i64_i64(int64_t a0, int64_t e0);
int64_t square_sum_below_power10_i32(int32_t k);
int64_t recurse_i64_i32_i32(int64_t suffix0, int32_t width, int32_t trailing_zeros);
int64_t correction_sum_i32(int32_t k);
int32_t main(void);
static const int64_t MOD = 1000000007;
/* Module statics */
static __int128* g_pow10 = ((__int128*)(NULL));
static int64_t* g_pow10_mod = ((int64_t*)(NULL));
static int32_t g_K = 0;
int64_t delta_mod_i128_i128(__int128 n0) {
if (n0 == 0) {
return 0;
}
int32_t digits[40] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
int32_t len = 0;
__int128 tmp = n0;
while (tmp > 0) {
digits[len] = ((int32_t)(FLOW_CHECKED_MOD((tmp), (10))));
len = (len + 1);
tmp = FLOW_CHECKED_DIV((tmp), (10));
}
int32_t i = 0;
int32_t j = (len - 1);
while (i < j) {
int32_t t = (((unsigned)(i) < 40) ? digits[i] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(i), 40), flow_fault_handler("array index out of bounds"), digits[0]));
digits[i] = (((unsigned)(j) < 40) ? digits[j] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(j), 40), flow_fault_handler("array index out of bounds"), digits[0]));
digits[j] = t;
i = (i + 1);
j = (j - 1);
}
int64_t n_mod = 0;
int32_t __flow_step_1 = 1;
for (int32_t d = 0; (0 <= len) ? d < len : d > len; d += (0 <= len) ? 1 : -1) {
n_mod = FLOW_CHECKED_MOD(((FLOW_CHECKED_MOD(((n_mod * 10)), (MOD)) + ((int64_t)((((unsigned)(d) < 40) ? digits[d] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(d), 40), flow_fault_handler("array index out of bounds"), digits[0])))))), (MOD));
}
int32_t ii = (len - 2);
while ((ii >= 0 && (((unsigned)(ii) < 40) ? digits[ii] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(ii), 40), flow_fault_handler("array index out of bounds"), digits[0])) >= (((unsigned)((ii + 1)) < 40) ? digits[(ii + 1)] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)((ii + 1)), 40), flow_fault_handler("array index out of bounds"), digits[0])))) {
ii = (ii - 1);
}
if (ii < 0) {
int64_t r = FLOW_CHECKED_MOD(((MOD - n_mod)), (MOD));
return r;
}
int32_t jj = (len - 1);
while ((((unsigned)(jj) < 40) ? digits[jj] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(jj), 40), flow_fault_handler("array index out of bounds"), digits[0])) <= (((unsigned)(ii) < 40) ? digits[ii] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(ii), 40), flow_fault_handler("array index out of bounds"), digits[0]))) {
jj = (jj - 1);
}
int32_t t = (((unsigned)(ii) < 40) ? digits[ii] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(ii), 40), flow_fault_handler("array index out of bounds"), digits[0]));
digits[ii] = (((unsigned)(jj) < 40) ? digits[jj] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(jj), 40), flow_fault_handler("array index out of bounds"), digits[0]));
digits[jj] = t;
int32_t lo = (ii + 1);
int32_t hi = (len - 1);
while (lo < hi) {
int32_t tt = (((unsigned)(lo) < 40) ? digits[lo] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(lo), 40), flow_fault_handler("array index out of bounds"), digits[0]));
digits[lo] = (((unsigned)(hi) < 40) ? digits[hi] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(hi), 40), flow_fault_handler("array index out of bounds"), digits[0]));
digits[hi] = tt;
lo = (lo + 1);
hi = (hi - 1);
}
int64_t np_mod = 0;
int32_t __flow_step_2 = 1;
for (int32_t d = 0; (0 <= len) ? d < len : d > len; d += (0 <= len) ? 1 : -1) {
np_mod = FLOW_CHECKED_MOD(((FLOW_CHECKED_MOD(((np_mod * 10)), (MOD)) + ((int64_t)((((unsigned)(d) < 40) ? digits[d] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(d), 40), flow_fault_handler("array index out of bounds"), digits[0])))))), (MOD));
}
int64_t dm = FLOW_CHECKED_MOD(((np_mod - n_mod)), (MOD));
if (dm < 0) {
return (dm + MOD);
}
return dm;
}
int64_t mod_pow_i64_i64(int64_t a0, int64_t e0) {
int64_t r = 1;
int64_t a = FLOW_CHECKED_MOD((a0), (MOD));
if (a < 0) {
a = (a + MOD);
}
int64_t e = e0;
while (e > 0) {
if ((e & 1) != 0) {
r = ((int64_t)(FLOW_CHECKED_MOD(((r * a)), (MOD))));
}
a = ((int64_t)(FLOW_CHECKED_MOD(((a * a)), (MOD))));
e = FLOW_CHECKED_SHR((e), (1));
}
return r;
}
int64_t square_sum_below_power10_i32(int32_t k) {
int64_t n = mod_pow_i64_i64(10, ((int64_t)(k)));
int64_t t1 = FLOW_CHECKED_MOD(((n - 1)), (MOD));
if (t1 < 0) {
t1 = (t1 + MOD);
}
int64_t t2 = FLOW_CHECKED_MOD((((2 * n) - 1)), (MOD));
if (t2 < 0) {
t2 = (t2 + MOD);
}
int64_t inv6 = mod_pow_i64_i64(6, (MOD - 2));
int64_t result = ((int64_t)(FLOW_CHECKED_MOD(((n * t1)), (MOD))));
result = ((int64_t)(FLOW_CHECKED_MOD(((result * t2)), (MOD))));
result = ((int64_t)(FLOW_CHECKED_MOD(((result * inv6)), (MOD))));
return result;
}
int64_t recurse_i64_i32_i32(int64_t suffix0, int32_t width, int32_t trailing_zeros) {
int32_t k = g_K;
if ((width + trailing_zeros) >= k) {
__int128 val = (((__int128)(suffix0)) * g_pow10[trailing_zeros]);
__int128 sq = (val * val);
return delta_mod_i128_i128(sq);
}
int64_t square_suffix = ((int64_t)(FLOW_CHECKED_MOD(((((__int128)(suffix0)) * ((__int128)(suffix0)))), (g_pow10[width]))));
int64_t left_digit = FLOW_CHECKED_DIV((square_suffix), (((int64_t)(g_pow10[(width - 1)]))));
__int128 next_place = g_pow10[width];
__int128 next_modulus = g_pow10[(width + 1)];
int32_t free_digits = (((k - width) - trailing_zeros) - 1);
int64_t completion_count = g_pow10_mod[free_digits];
int64_t total = 0;
int32_t __flow_step_3 = 1;
for (int32_t digit = 0; (0 <= 10) ? digit < 10 : digit > 10; digit += (0 <= 10) ? 1 : -1) {
__int128 next_suffix = ((((__int128)(digit)) * next_place) + ((__int128)(suffix0)));
__int128 next_sq = (next_suffix * next_suffix);
int64_t new_digit = ((int64_t)(FLOW_CHECKED_DIV((FLOW_CHECKED_MOD((next_sq), (next_modulus))), (next_place))));
if (new_digit >= left_digit) {
total = FLOW_CHECKED_MOD(((total + recurse_i64_i32_i32(((int64_t)(next_suffix)), (width + 1), trailing_zeros))), (MOD));
} else {
__int128 next_sq_mod = FLOW_CHECKED_MOD((next_sq), (next_modulus));
__int128 visible = (next_sq_mod * g_pow10[(2 * trailing_zeros)]);
__int128 representative = (g_pow10[((width + 1) + (2 * trailing_zeros))] + visible);
int64_t representative_delta = delta_mod_i128_i128(representative);
if (next_sq < next_modulus) {
total = FLOW_CHECKED_MOD(((total + delta_mod_i128_i128(visible))), (MOD));
int64_t term = FLOW_CHECKED_MOD(((completion_count - 1)), (MOD));
if (term < 0) {
term = (term + MOD);
}
term = ((int64_t)(FLOW_CHECKED_MOD(((term * representative_delta)), (MOD))));
total = FLOW_CHECKED_MOD(((total + term)), (MOD));
} else {
int64_t term = ((int64_t)(FLOW_CHECKED_MOD(((completion_count * representative_delta)), (MOD))));
total = FLOW_CHECKED_MOD(((total + term)), (MOD));
}
}
}
return total;
}
int64_t correction_sum_i32(int32_t k) {
g_K = k;
g_pow10[0] = 1;
int32_t __flow_step_4 = 1;
for (int32_t i = 1; (1 <= ((2 * k) + 3)) ? i < ((2 * k) + 3) : i > ((2 * k) + 3); i += (1 <= ((2 * k) + 3)) ? 1 : -1) {
g_pow10[i] = (g_pow10[(i - 1)] * 10);
}
g_pow10_mod[0] = 1;
int32_t __flow_step_5 = 1;
for (int32_t i = 1; (1 <= (k + 1)) ? i < (k + 1) : i > (k + 1); i += (1 <= (k + 1)) ? 1 : -1) {
g_pow10_mod[i] = FLOW_CHECKED_MOD(((g_pow10_mod[(i - 1)] * 10)), (MOD));
}
int64_t total = 0;
int32_t __flow_step_6 = 1;
for (int32_t trailing_zeros = 0; (0 <= k) ? trailing_zeros < k : trailing_zeros > k; trailing_zeros += (0 <= k) ? 1 : -1) {
int32_t __flow_step_7 = 1;
for (int32_t last_digit = 1; (1 <= 10) ? last_digit < 10 : last_digit > 10; last_digit += (1 <= 10) ? 1 : -1) {
total = FLOW_CHECKED_MOD(((total + recurse_i64_i32_i32(((int64_t)(last_digit)), 1, trailing_zeros))), (MOD));
}
}
return total;
}
int32_t main(void) {
g_pow10 = ((__int128*)(malloc((40 * 16))));
g_pow10_mod = ((int64_t*)(malloc((20 * 8))));
int32_t k = 16;
int64_t result = FLOW_CHECKED_MOD(((square_sum_below_power10_i32(k) + correction_sum_i32(k))), (MOD));
printf("%lld\n", result);
return 0;
}