Inverse Digit Sum II: S(10^4) = sum f(n^3, n^4) mod 1e9+7, where f(s, m) is the m-th positive integer with digit sum s. Numbers with digit sum s and d digits are all-9 strings with a small total "deficiency" t = 9d - s spread over the digits (first digit loses at most 8). Choose the digit count by comparing m against per-length counts, then unrank digit by digit: at each position try the largest deficiency first, counting completions with capped-at-1e30 i128 binomials (any capped count exceeds every possible m, so comparisons stay exact). Runs of unchanged 9s are jumped in closed form. The value mod p is 10^d - 1 minus the few lowered digits.
# Project Euler 685
# Inverse Digit Sum II: S(10^4) = sum f(n^3, n^4) mod 1e9+7, where f(s, m)
# is the m-th positive integer with digit sum s.
#
# Numbers with digit sum s and d digits are all-9 strings with a small
# total "deficiency" t = 9d - s spread over the digits (first digit loses
# at most 8). Choose the digit count by comparing m against per-length
# counts, then unrank digit by digit: at each position try the largest
# deficiency first, counting completions with capped-at-1e30 i128
# binomials (any capped count exceeds every possible m, so comparisons
# stay exact). Runs of unchanged 9s are jumped in closed form. The value
# mod p is 10^d - 1 minus the few lowered digits.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
const MOD: i64 = 1000000007
const KMAX: i64 = 10000
function powmod(base: i64, exp: i64) -> i64 {
let mut result: i64 = 1
let mut b: i64 = base % MOD
let mut e: i64 = exp
while e > 0 {
if (e & 1) == 1 {
result = result * b % MOD
}
b = b * b % MOD
e = e >> 1
}
return result
}
function cap128() -> i128 {
let a: i128 = 1000000000000000 as i128
return a * a # 1e30
}
# C(n, k) exact in i128, saturating at 1e30
function binom_capped(n: i64, k0: i64) -> i128 {
if k0 < 0 || n < 0 || k0 > n {
return 0 as i128
}
let mut k: i64 = k0
if n - k < k {
k = n - k
}
let CAP: i128 = cap128()
let mut c: i128 = 1 as i128
let mut i: i64 = 1
while i <= k {
let fac: i128 = (n - k + i) as i128
if c > CAP / fac {
return CAP
}
c = c * fac / (i as i128)
if c >= CAP {
return CAP
}
i = i + 1
}
return c
}
# number of ways to write tau as rho digits each in [0,9]; capped at 1e30
function Wcnt(rho: i64, tau: i64) -> i128 {
if tau < 0 || tau > 9 * rho {
return 0 as i128
}
if rho == 0 {
if tau == 0 {
return 1 as i128
}
return 0 as i128
}
let CAP: i128 = cap128()
let mut tot: i128 = 0 as i128
let mut sign: i64 = 1
let mut k: i64 = 0
let mut t2: i64 = tau
while t2 >= 0 {
let term: i128 = binom_capped(rho, k) * 1 as i128 * binom_capped(t2 + rho - 1, t2)
# (binom_capped(rho,k) is tiny for k<=3; product cannot overflow the cap regime)
if term >= CAP {
# capped counts always exceed any m; sign is + here (k=0 dominates)
return CAP
}
if sign == 1 {
tot = tot + term
} else {
tot = tot - term
}
sign = 0 - sign
k = k + 1
t2 = t2 - 10
}
if tot < 0 as i128 {
tot = 0 as i128
}
if tot > CAP {
tot = CAP
}
return tot
}
# d-digit numbers with deficiency t (first digit >= 1)
function Wfull(d: i64, t: i64) -> i128 {
let CAP: i128 = cap128()
let mut tot: i128 = 0 as i128
let mut b1: i64 = 0
while b1 <= 8 && b1 <= t {
tot = tot + Wcnt(d - 1, t - b1)
if tot > CAP {
tot = CAP
}
b1 = b1 + 1
}
return tot
}
# f(s, m) mod MOD
function f_mod(s: i64, m0: i64) -> i64 {
let mut d: i64 = (s + 8) / 9
let mut t: i64 = 9 * d - s
let mut m: i64 = m0
let mut choosing: i64 = 1
while choosing == 1 {
let tot: i128 = Wfull(d, t)
if (m as i128) > tot {
m = m - (tot as i64)
d = d + 1
t = t + 9
} else {
choosing = 0
}
}
# unrank: sparse deviations (position, lowered-by)
let dpos: ptr<i64> = calloc(64, 8)
let dval: ptr<i64> = calloc(64, 8)
let mut nd: i64 = 0
let mut i: i64 = 1
let mut trem: i64 = t
while trem > 0 {
let r: i64 = d - i + 1
let mut cap: i64 = 9
if i == 1 {
cap = 8
}
let mut placed: i64 = 0
let mut b: i64 = trem
if b > cap {
b = cap
}
while b >= 1 && placed == 0 {
let w: i128 = Wcnt(r - 1, trem - b)
if (m as i128) > w {
m = m - (w as i64)
b = b - 1
} else {
dpos[nd] = i
dval[nd] = b
nd = nd + 1
trem = trem - b
i = i + 1
placed = 1
}
}
if placed == 0 {
# this position keeps digit 9; jump the run
if trem == 1 {
# the m-th remaining vector has its single 1 at position i+m
dpos[nd] = i + m
dval[nd] = 1
nd = nd + 1
trem = 0
} else {
if trem == 2 {
# cum(q) = sum_{j=1..q} (r-j) = q*r - q(q+1)/2 in i128
let mut lo: i64 = 1
let mut hi: i64 = r - 1
while lo < hi {
let mid: i64 = (lo + hi) / 2
let cum: i128 = (mid as i128) * (r as i128)
- (mid as i128) * ((mid + 1) as i128) / (2 as i128)
if cum >= (m as i128) {
hi = mid
} else {
lo = mid + 1
}
}
let q: i64 = lo
let used: i128 = ((q - 1) as i128) * (r as i128)
- ((q - 1) as i128) * (q as i128) / (2 as i128)
m = m - (used as i64)
i = i + q
} else {
i = i + 1
}
}
}
}
# value = (10^d - 1) - sum b * 10^(d - pos) (mod MOD)
let mut v: i64 = (powmod(10, d) - 1 + MOD) % MOD
let mut u: i64 = 0
while u < nd {
v = (v - dval[u] * powmod(10, d - dpos[u]) % MOD + MOD) % MOD
u = u + 1
}
free(dpos)
free(dval)
return v
}
function main() -> i32 {
let mut S: i64 = 0
let mut n: i64 = 1
while n <= KMAX {
let s: i64 = n * n * n
let m: i64 = s * n
S = (S + f_mod(s, m)) % MOD
n = n + 1
}
printf("%lld\n", S)
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 powmod_i64_i64(int64_t base, int64_t exp);
__int128 cap128(void);
__int128 binom_capped_i64_i64(int64_t n, int64_t k0);
__int128 Wcnt_i64_i64(int64_t rho, int64_t tau);
__int128 Wfull_i64_i64(int64_t d, int64_t t);
int64_t f_mod_i64_i64(int64_t s, int64_t m0);
int32_t main(void);
static const int64_t MOD = 1000000007;
static const int64_t KMAX = 10000;
int64_t powmod_i64_i64(int64_t base, int64_t exp) {
int64_t result = 1;
int64_t b = FLOW_CHECKED_MOD((base), (MOD));
int64_t e = exp;
while (e > 0) {
if ((e & 1) == 1) {
result = FLOW_CHECKED_MOD(((result * b)), (MOD));
}
b = FLOW_CHECKED_MOD(((b * b)), (MOD));
e = FLOW_CHECKED_SHR((e), (1));
}
return result;
}
__int128 cap128(void) {
__int128 a = ((__int128)(1000000000000000));
return (a * a);
}
__int128 binom_capped_i64_i64(int64_t n, int64_t k0) {
if (((k0 < 0 || n < 0) || k0 > n)) {
return ((__int128)(0));
}
int64_t k = k0;
if ((n - k) < k) {
k = (n - k);
}
__int128 CAP = cap128();
__int128 c = ((__int128)(1));
int64_t i = 1;
while (i <= k) {
__int128 fac = ((__int128)(((n - k) + i)));
if (c > FLOW_CHECKED_DIV((CAP), (fac))) {
return CAP;
}
c = FLOW_CHECKED_DIV(((c * fac)), (((__int128)(i))));
if (c >= CAP) {
return CAP;
}
i = (i + 1);
}
return c;
}
__int128 Wcnt_i64_i64(int64_t rho, int64_t tau) {
if ((tau < 0 || tau > (9 * rho))) {
return ((__int128)(0));
}
if (rho == 0) {
if (tau == 0) {
return ((__int128)(1));
}
return ((__int128)(0));
}
__int128 CAP = cap128();
__int128 tot = ((__int128)(0));
int64_t sign = 1;
int64_t k = 0;
int64_t t2 = tau;
while (t2 >= 0) {
__int128 term = ((binom_capped_i64_i64(rho, k) * ((__int128)(1))) * binom_capped_i64_i64(((t2 + rho) - 1), t2));
if (term >= CAP) {
return CAP;
}
if (sign == 1) {
tot = (tot + term);
} else {
tot = (tot - term);
}
sign = (0 - sign);
k = (k + 1);
t2 = (t2 - 10);
}
if (tot < ((__int128)(0))) {
tot = ((__int128)(0));
}
if (tot > CAP) {
tot = CAP;
}
return tot;
}
__int128 Wfull_i64_i64(int64_t d, int64_t t) {
__int128 CAP = cap128();
__int128 tot = ((__int128)(0));
int64_t b1 = 0;
while ((b1 <= 8 && b1 <= t)) {
tot = (tot + Wcnt_i64_i64((d - 1), (t - b1)));
if (tot > CAP) {
tot = CAP;
}
b1 = (b1 + 1);
}
return tot;
}
int64_t f_mod_i64_i64(int64_t s, int64_t m0) {
int64_t d = FLOW_CHECKED_DIV(((s + 8)), (9));
int64_t t = ((9 * d) - s);
int64_t m = m0;
int64_t choosing = 1;
while (choosing == 1) {
__int128 tot = Wfull_i64_i64(d, t);
if (((__int128)(m)) > tot) {
m = (m - ((int64_t)(tot)));
d = (d + 1);
t = (t + 9);
} else {
choosing = 0;
}
}
int64_t* dpos = (int64_t*)(calloc(64, 8));
int64_t* dval = (int64_t*)(calloc(64, 8));
int64_t nd = 0;
int64_t i = 1;
int64_t trem = t;
while (trem > 0) {
int64_t r = ((d - i) + 1);
int64_t cap = 9;
if (i == 1) {
cap = 8;
}
int64_t placed = 0;
int64_t b = trem;
if (b > cap) {
b = cap;
}
while ((b >= 1 && placed == 0)) {
__int128 w = Wcnt_i64_i64((r - 1), (trem - b));
if (((__int128)(m)) > w) {
m = (m - ((int64_t)(w)));
b = (b - 1);
} else {
dpos[nd] = i;
dval[nd] = b;
nd = (nd + 1);
trem = (trem - b);
i = (i + 1);
placed = 1;
}
}
if (placed == 0) {
if (trem == 1) {
dpos[nd] = (i + m);
dval[nd] = 1;
nd = (nd + 1);
trem = 0;
} else {
if (trem == 2) {
int64_t lo = 1;
int64_t hi = (r - 1);
while (lo < hi) {
int64_t mid = FLOW_CHECKED_DIV(((lo + hi)), (2));
__int128 cum = ((((__int128)(mid)) * ((__int128)(r))) - FLOW_CHECKED_DIV(((((__int128)(mid)) * ((__int128)((mid + 1))))), (((__int128)(2)))));
if (cum >= ((__int128)(m))) {
hi = mid;
} else {
lo = (mid + 1);
}
}
int64_t q = lo;
__int128 used = ((((__int128)((q - 1))) * ((__int128)(r))) - FLOW_CHECKED_DIV(((((__int128)((q - 1))) * ((__int128)(q)))), (((__int128)(2)))));
m = (m - ((int64_t)(used)));
i = (i + q);
} else {
i = (i + 1);
}
}
}
}
int64_t v = FLOW_CHECKED_MOD((((powmod_i64_i64(10, d) - 1) + MOD)), (MOD));
int64_t u = 0;
while (u < nd) {
v = FLOW_CHECKED_MOD((((v - FLOW_CHECKED_MOD(((dval[u] * powmod_i64_i64(10, (d - dpos[u])))), (MOD))) + MOD)), (MOD));
u = (u + 1);
}
free(dpos);
free(dval);
return v;
}
int32_t main(void) {
int64_t S = 0;
int64_t n = 1;
while (n <= KMAX) {
int64_t s = ((n * n) * n);
int64_t m = (s * n);
S = FLOW_CHECKED_MOD(((S + f_mod_i64_i64(s, m))), (MOD));
n = (n + 1);
}
printf("%lld\n", S);
return 0;
}