Best Approximations by Quadratic Integers — i128 fixed-point (96 fractional bits). Machin's formula for pi, Newton's method for sqrt, continued fraction best approximation.
# Project Euler 591
# Best Approximations by Quadratic Integers — i128 fixed-point (96 fractional bits).
# Machin's formula for pi, Newton's method for sqrt, continued fraction best approximation.
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
function realloc(p: ptr<void>, n: i64) -> ptr<void>
function printf(fmt: ptr<i8>, ...) -> i32
function sqrt(x: f64) -> f64
function floor(x: f64) -> f64
function round(x: f64) -> f64
}
function one() -> i128 {
return (1 as i128) << 96
}
function arctan_inv_scaled(x: i32) -> i128 {
let mut sum: i128 = 0
let mut xpow: i128 = x as i128
let x2: i128 = (x as i128) * (x as i128)
let mut k: i32 = 0
while k < 60 {
let denom: i128 = ((2 * k + 1) as i128) * xpow
if denom > one() { break }
let term: i128 = one() / denom
if k % 2 == 0 {
sum = sum + term
} else {
sum = sum - term
}
xpow = xpow * x2
k = k + 1
}
return sum
}
function compute_pi_fp() -> i128 {
return 16 * arctan_inv_scaled(5) - 4 * arctan_inv_scaled(239)
}
function floor_div(a: i128, b: i128) -> i128 {
let q: i128 = a / b
let r: i128 = a % b
if r != 0 {
if (a < 0 && b >= 0) || (a >= 0 && b < 0) {
return q - 1
}
}
return q
}
function isqrt_scaled(d: i32) -> i128 {
let target: i128 = (d as i128) << 120
let shift60: f64 = ((1 as i64) << 60) as f64
let mut y: i128 = (sqrt(d as f64) * shift60) as i128
let mut i: i32 = 0
while i < 50 {
let q: i128 = target / y
let y_new: i128 = (y + q) / 2
if y_new == y { break }
y = y_new
i = i + 1
}
let mut z: i128 = y << 36
let mut two_sqrt_d: i32 = round(2.0 * sqrt(d as f64)) as i32
if two_sqrt_d == 0 { two_sqrt_d = 1 }
let mut iter: i32 = 0
while iter < 40 {
let hi: i128 = z >> 48
let lo: i128 = z & (((1 as i128) << 48) - 1)
let z2: i128 = hi * hi + ((2 * hi * lo) >> 48) + ((lo * lo) >> 96)
let diff: i128 = z2 - (d as i128) * one()
if diff == 0 { break }
let mut corr: i128 = 0
if diff >= 0 - ((1 as i128) << 30) && diff <= ((1 as i128) << 30) {
corr = floor_div(diff * one(), (2 as i128) * z)
if corr == 0 {
if diff > 0 { corr = 1 } else { corr = 0 - 1 }
}
} else {
corr = floor_div(diff, two_sqrt_d as i128)
}
if corr == 0 { break }
z = z - corr
iter = iter + 1
}
return z
}
function frac_mul(alpha_fp: i128, n: i64) -> i128 {
let mask64: i128 = ((1 as i128) << 64) - 1
let lo: i128 = alpha_fp & mask64
let hi: i128 = (alpha_fp >> 64) & mask64
let lo_n: i128 = lo * (n as i128)
let hi_n: i128 = hi * (n as i128)
let hi_n_mod: i128 = hi_n & (((1 as i128) << 32) - 1)
let hi_part: i128 = hi_n_mod << 64
return (lo_n + hi_part) & (one() - 1)
}
function compute_error_fp(sqrt_D_fp: i128, pi_fp: i128, b: i64) -> i128 {
let mask64: i128 = ((1 as i128) << 64) - 1
let lo_sqrt: i128 = sqrt_D_fp & mask64
let hi_sqrt: i128 = (sqrt_D_fp >> 64) & mask64
let lo_b: i128 = lo_sqrt * (b as i128)
let hi_b: i128 = hi_sqrt * (b as i128)
let hi_b_mod: i128 = hi_b & (((1 as i128) << 32) - 1)
let hi_part: i128 = hi_b_mod << 64
let val: i128 = lo_b - pi_fp + hi_part
let mut frac_val: i128 = val % one()
if frac_val < 0 { frac_val = frac_val + one() }
if frac_val <= one() / 2 { return frac_val }
return one() - frac_val
}
function is_square(n: i32) -> bool {
let r: i32 = floor(sqrt(n as f64) + 0.5) as i32
return r * r == n
}
function sqrt_cf_period(d: i32, a0_out: ptr<i32>, period: ptr<i32>, plen_out: ptr<i32>) -> i32 {
let a0: i32 = floor(sqrt(d as f64)) as i32
a0_out[0] = a0
let mut m: i64 = 0
let mut dd: i64 = 1
let mut a: i64 = a0 as i64
let mut n: i32 = 0
while true {
m = dd * a - m
dd = ((d as i64) - m * m) / dd
a = (a0 as i64 + m) / dd
period[n] = a as i32
n = n + 1
if a == 2 * (a0 as i64) { break }
if n > 10000 { return -1 }
}
plen_out[0] = n
return 0
}
function ceil_div_i128(a0: i128, b0: i128) -> i128 {
let mut a: i128 = a0
let mut b: i128 = b0
if b < 0 {
a = 0 - a
b = 0 - b
}
if a >= 0 { return (a + b - 1) / b }
return 0 - ((0 - a) / b)
}
function best_b_positive(alpha_fp: i128, beta_fp: i128, B: i64, period: ptr<i32>, plen: i32) -> i64 {
if B <= 0 { return 0 }
let mut cap: i32 = 128
let mut a: ptr<i32> = calloc(cap as i64, 4) as ptr<i32>
let mut q: ptr<i64> = calloc(cap as i64, 8) as ptr<i64>
let mut delta: ptr<i128> = calloc(cap as i64, 16) as ptr<i128>
a[0] = 0
q[0] = 1
let mut q_minus1: i64 = 0
delta[0] = alpha_fp
let delta_minus1: i128 = one()
let mut k: i32 = 1
let mut extra: i32 = 6
let mut len: i32 = 1
while true {
if k + 2 >= cap {
let ncap: i32 = cap * 2
a = realloc(a as ptr<void>, (ncap as i64) * 4) as ptr<i32>
q = realloc(q as ptr<void>, (ncap as i64) * 8) as ptr<i64>
delta = realloc(delta as ptr<void>, (ncap as i64) * 16) as ptr<i128>
cap = ncap
}
let ak: i32 = period[(k - 1) % plen]
a[k] = ak
let qk: i64 = (ak as i64) * q[k - 1] + q_minus1
q_minus1 = q[k - 1]
q[k] = qk
if k == 1 {
delta[k] = 0 - (ak as i128) * delta[0] + delta_minus1
} else {
delta[k] = 0 - (ak as i128) * delta[k - 1] + delta[k - 2]
}
len = k + 1
if qk > B {
extra = extra - 1
if extra <= 0 { break }
}
k = k + 1
if k > 500 { break }
}
let max_i: i32 = len - 1
let b_digits: ptr<i32> = calloc((max_i + 1) as i64, 4) as ptr<i32>
let mut beta_rem: i128 = beta_fp
let mut i: i32 = 1
while i <= max_i {
let d: i128 = delta[i - 1]
if d == 0 {
b_digits[i] = 0
} else {
let mut bi: i128 = ceil_div_i128(beta_rem, d)
if bi > (a[i] as i128) { bi = a[i] as i128 }
if bi < 0 { bi = 0 }
b_digits[i] = bi as i32
beta_rem = bi * d - beta_rem
}
i = i + 1
}
let prefix: ptr<i64> = calloc((max_i + 1) as i64, 8) as ptr<i64>
let mut s: i64 = 0
i = 1
while i <= max_i {
s = s + (b_digits[i] as i64) * q[i - 1]
prefix[i] = s
i = i + 1
}
let cand_r: ptr<i64> = calloc(200000, 8) as ptr<i64>
let cand_l: ptr<i64> = calloc(200000, 8) as ptr<i64>
let mut nr: i32 = 0
let mut nl: i32 = 0
cand_r[nr] = 0
nr = nr + 1
cand_l[nl] = 0
nl = nl + 1
let mut k2: i32 = 1
while 2 * k2 < max_i + 1 {
let ie: i32 = 2 * k2
let io: i32 = 2 * k2 - 1
if ie > max_i { break }
let P: i64 = prefix[io]
let st: i64 = q[io]
let mut j: i32 = 0
while j < b_digits[ie] {
let n: i64 = P + (j as i64) * st
if 0 <= n && n <= B {
cand_r[nr] = n
nr = nr + 1
}
j = j + 1
}
k2 = k2 + 1
}
k2 = 0
while 2 * k2 + 1 <= max_i {
let idx: i32 = 2 * k2
let inn: i32 = idx + 1
if inn > max_i { break }
let P: i64 = prefix[idx]
let st: i64 = q[idx]
let mut j: i32 = 0
while j < b_digits[inn] {
let n: i64 = P + (j as i64) * st
if 0 <= n && n <= B {
cand_l[nl] = n
nl = nl + 1
}
j = j + 1
}
k2 = k2 + 1
}
let mut best_r_gap: i128 = one()
let mut best_l_gap: i128 = one()
let mut best_r_n: i64 = 0
let mut best_l_n: i64 = 0
i = 0
while i < nr {
let n: i64 = cand_r[i]
let x: i128 = frac_mul(alpha_fp, n)
let mut gap: i128 = x - beta_fp
if gap < 0 { gap = gap + one() }
if gap < best_r_gap {
best_r_gap = gap
best_r_n = n
}
i = i + 1
}
i = 0
while i < nl {
let n: i64 = cand_l[i]
let x: i128 = frac_mul(alpha_fp, n)
let mut gap: i128 = beta_fp - x
if gap < 0 { gap = gap + one() }
if gap < best_l_gap {
best_l_gap = gap
best_l_n = n
}
i = i + 1
}
let ans: i64 = 0
if best_r_gap < best_l_gap {
ans = best_r_n
} else {
ans = best_l_n
}
free(delta as ptr<void>)
free(a as ptr<void>)
free(q as ptr<void>)
free(b_digits as ptr<void>)
free(prefix as ptr<void>)
free(cand_r as ptr<void>)
free(cand_l as ptr<void>)
return ans
}
function bqa_pi_d(d: i32, n: i64, pi_fp: i128, beta_fp: i128) -> i64 {
let a0_ptr: ptr<i32> = calloc(1, 4) as ptr<i32>
let period: ptr<i32> = calloc(4096, 4) as ptr<i32>
let plen_ptr: ptr<i32> = calloc(1, 4) as ptr<i32>
sqrt_cf_period(d, a0_ptr, period, plen_ptr)
let a0: i32 = a0_ptr[0]
let plen: i32 = plen_ptr[0]
let sqrt_D_fp: i128 = isqrt_scaled(d)
let alpha_fp: i128 = sqrt_D_fp - ((a0 as i128) << 96)
let sqrt_d: f64 = sqrt(d as f64)
let pi_d: f64 = 3.14159265358979323846
let mut Bpos: i64 = floor((pi_d + (n as f64)) / sqrt_d) as i64
if Bpos < 0 { Bpos = 0 }
if Bpos > n { Bpos = n }
let mut Bneg: i64 = floor(((n as f64) - pi_d) / sqrt_d) as i64
if Bneg < 0 { Bneg = 0 }
if Bneg > n { Bneg = n }
let b_pos: i64 = best_b_positive(alpha_fp, beta_fp, Bpos, period, plen)
let one_m_beta: i128 = one() - beta_fp
let t: i64 = best_b_positive(alpha_fp, one_m_beta, Bneg, period, plen)
let b_neg: i64 = 0 - t
let err1_fp: i128 = compute_error_fp(sqrt_D_fp, pi_fp, b_pos)
let err2_fp: i128 = compute_error_fp(sqrt_D_fp, pi_fp, b_neg)
let mut b_chosen: i64 = b_pos
if err2_fp < err1_fp { b_chosen = b_neg }
let a_real: f64 = pi_d - sqrt_d * (b_chosen as f64)
let mut a_out: i64 = round(a_real) as i64
if a_out > n { a_out = n }
if a_out < 0 - n { a_out = 0 - n }
free(a0_ptr as ptr<void>)
free(period as ptr<void>)
free(plen_ptr as ptr<void>)
return a_out
}
function main() -> i32 {
let pi_fp: i128 = compute_pi_fp()
let beta_fp: i128 = pi_fp - ((3 as i128) << 96)
let mut total: i64 = 0
let n: i64 = 10000000000000
let mut d: i32 = 2
while d < 100 {
if !is_square(d) {
let a: i64 = bqa_pi_d(d, n, pi_fp, beta_fp)
if a >= 0 {
total = total + a
} else {
total = total + (0 - a)
}
}
d = d + 1
}
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; }
__int128 one(void);
__int128 arctan_inv_scaled_i32(int32_t x);
__int128 compute_pi_fp(void);
__int128 floor_div_i128_i128(__int128 a, __int128 b);
__int128 isqrt_scaled_i32(int32_t d);
__int128 frac_mul_i128_i64(__int128 alpha_fp, int64_t n);
__int128 compute_error_fp_i128_i128_i64(__int128 sqrt_D_fp, __int128 pi_fp, int64_t b);
bool is_square_i32(int32_t n);
int32_t sqrt_cf_period_i32_ptr_i32_ptr_i32_ptr_i32(int32_t d, int32_t* a0_out, int32_t* period, int32_t* plen_out);
__int128 ceil_div_i128_i128_i128(__int128 a0, __int128 b0);
int64_t best_b_positive_i128_i128_i64_ptr_i32_i32(__int128 alpha_fp, __int128 beta_fp, int64_t B, int32_t* period, int32_t plen);
int64_t bqa_pi_d_i32_i64_i128_i128(int32_t d, int64_t n, __int128 pi_fp, __int128 beta_fp);
int32_t main(void);
__int128 one(void) {
return FLOW_CHECKED_SHL((((__int128)(1))), (96));
}
__int128 arctan_inv_scaled_i32(int32_t x) {
__int128 sum = 0;
__int128 xpow = ((__int128)(x));
__int128 x2 = (((__int128)(x)) * ((__int128)(x)));
int32_t k = 0;
while (k < 60) {
__int128 denom = (((__int128)(((2 * k) + 1))) * xpow);
if (denom > one()) {
break;
}
__int128 term = FLOW_CHECKED_DIV((one()), (denom));
if (FLOW_CHECKED_MOD((k), (2)) == 0) {
sum = (sum + term);
} else {
sum = (sum - term);
}
xpow = (xpow * x2);
k = (k + 1);
}
return sum;
}
__int128 compute_pi_fp(void) {
return ((16 * arctan_inv_scaled_i32(5)) - (4 * arctan_inv_scaled_i32(239)));
}
__int128 floor_div_i128_i128(__int128 a, __int128 b) {
__int128 q = FLOW_CHECKED_DIV((a), (b));
__int128 r = FLOW_CHECKED_MOD((a), (b));
if (r != 0) {
if (((a < 0 && b >= 0) || (a >= 0 && b < 0))) {
return (q - 1);
}
}
return q;
}
__int128 isqrt_scaled_i32(int32_t d) {
__int128 target = FLOW_CHECKED_SHL((((__int128)(d))), (120));
double shift60 = ((double)(FLOW_CHECKED_SHL((((int64_t)(1))), (60))));
__int128 y = ((__int128)((sqrt(((double)(d))) * shift60)));
int32_t i = 0;
while (i < 50) {
__int128 q = FLOW_CHECKED_DIV((target), (y));
__int128 y_new = FLOW_CHECKED_DIV(((y + q)), (2));
if (y_new == y) {
break;
}
y = y_new;
i = (i + 1);
}
__int128 z = FLOW_CHECKED_SHL((y), (36));
int32_t two_sqrt_d = ((int32_t)(round((2.0 * sqrt(((double)(d)))))));
if (two_sqrt_d == 0) {
two_sqrt_d = 1;
}
int32_t iter = 0;
while (iter < 40) {
__int128 hi = FLOW_CHECKED_SHR((z), (48));
__int128 lo = (z & (FLOW_CHECKED_SHL((((__int128)(1))), (48)) - 1));
__int128 z2 = (((hi * hi) + FLOW_CHECKED_SHR((((2 * hi) * lo)), (48))) + FLOW_CHECKED_SHR(((lo * lo)), (96)));
__int128 diff = (z2 - (((__int128)(d)) * one()));
if (diff == 0) {
break;
}
__int128 corr = 0;
if ((diff >= (0 - FLOW_CHECKED_SHL((((__int128)(1))), (30))) && diff <= FLOW_CHECKED_SHL((((__int128)(1))), (30)))) {
corr = floor_div_i128_i128((diff * one()), (((__int128)(2)) * z));
if (corr == 0) {
if (diff > 0) {
corr = 1;
} else {
corr = (0 - 1);
}
}
} else {
corr = floor_div_i128_i128(diff, ((__int128)(two_sqrt_d)));
}
if (corr == 0) {
break;
}
z = (z - corr);
iter = (iter + 1);
}
return z;
}
__int128 frac_mul_i128_i64(__int128 alpha_fp, int64_t n) {
__int128 mask64 = (FLOW_CHECKED_SHL((((__int128)(1))), (64)) - 1);
__int128 lo = (alpha_fp & mask64);
__int128 hi = (FLOW_CHECKED_SHR((alpha_fp), (64)) & mask64);
__int128 lo_n = (lo * ((__int128)(n)));
__int128 hi_n = (hi * ((__int128)(n)));
__int128 hi_n_mod = (hi_n & (FLOW_CHECKED_SHL((((__int128)(1))), (32)) - 1));
__int128 hi_part = FLOW_CHECKED_SHL((hi_n_mod), (64));
return ((lo_n + hi_part) & (one() - 1));
}
__int128 compute_error_fp_i128_i128_i64(__int128 sqrt_D_fp, __int128 pi_fp, int64_t b) {
__int128 mask64 = (FLOW_CHECKED_SHL((((__int128)(1))), (64)) - 1);
__int128 lo_sqrt = (sqrt_D_fp & mask64);
__int128 hi_sqrt = (FLOW_CHECKED_SHR((sqrt_D_fp), (64)) & mask64);
__int128 lo_b = (lo_sqrt * ((__int128)(b)));
__int128 hi_b = (hi_sqrt * ((__int128)(b)));
__int128 hi_b_mod = (hi_b & (FLOW_CHECKED_SHL((((__int128)(1))), (32)) - 1));
__int128 hi_part = FLOW_CHECKED_SHL((hi_b_mod), (64));
__int128 val = ((lo_b - pi_fp) + hi_part);
__int128 frac_val = FLOW_CHECKED_MOD((val), (one()));
if (frac_val < 0) {
frac_val = (frac_val + one());
}
if (frac_val <= FLOW_CHECKED_DIV((one()), (2))) {
return frac_val;
}
return (one() - frac_val);
}
bool is_square_i32(int32_t n) {
int32_t r = ((int32_t)(floor((sqrt(((double)(n))) + 0.5))));
return (r * r) == n;
}
int32_t sqrt_cf_period_i32_ptr_i32_ptr_i32_ptr_i32(int32_t d, int32_t* a0_out, int32_t* period, int32_t* plen_out) {
int32_t a0 = ((int32_t)(floor(sqrt(((double)(d))))));
a0_out[0] = a0;
int64_t m = 0;
int64_t dd = 1;
int64_t a = ((int64_t)(a0));
int32_t n = 0;
while (1) {
m = ((dd * a) - m);
dd = FLOW_CHECKED_DIV(((((int64_t)(d)) - (m * m))), (dd));
a = FLOW_CHECKED_DIV(((((int64_t)(a0)) + m)), (dd));
period[n] = ((int32_t)(a));
n = (n + 1);
if (a == (2 * ((int64_t)(a0)))) {
break;
}
if (n > 10000) {
return (-1);
}
}
plen_out[0] = n;
return 0;
}
__int128 ceil_div_i128_i128_i128(__int128 a0, __int128 b0) {
__int128 a = a0;
__int128 b = b0;
if (b < 0) {
a = (0 - a);
b = (0 - b);
}
if (a >= 0) {
return FLOW_CHECKED_DIV((((a + b) - 1)), (b));
}
return (0 - FLOW_CHECKED_DIV(((0 - a)), (b)));
}
int64_t best_b_positive_i128_i128_i64_ptr_i32_i32(__int128 alpha_fp, __int128 beta_fp, int64_t B, int32_t* period, int32_t plen) {
if (B <= 0) {
return 0;
}
int32_t cap = 128;
int32_t* a = (int32_t*)(((int32_t*)(calloc(((int64_t)(cap)), 4))));
int64_t* q = (int64_t*)(((int64_t*)(calloc(((int64_t)(cap)), 8))));
__int128* delta = (__int128*)(((__int128*)(calloc(((int64_t)(cap)), 16))));
a[0] = 0;
q[0] = 1;
int64_t q_minus1 = 0;
delta[0] = alpha_fp;
__int128 delta_minus1 = one();
int32_t k = 1;
int32_t extra = 6;
int32_t len = 1;
while (1) {
if ((k + 2) >= cap) {
int32_t ncap = (cap * 2);
a = ((int32_t*)(realloc(((void*)(a)), (((int64_t)(ncap)) * 4))));
q = ((int64_t*)(realloc(((void*)(q)), (((int64_t)(ncap)) * 8))));
delta = ((__int128*)(realloc(((void*)(delta)), (((int64_t)(ncap)) * 16))));
cap = ncap;
}
int32_t ak = period[FLOW_CHECKED_MOD(((k - 1)), (plen))];
a[k] = ak;
int64_t qk = ((((int64_t)(ak)) * q[(k - 1)]) + q_minus1);
q_minus1 = q[(k - 1)];
q[k] = qk;
if (k == 1) {
delta[k] = ((0 - (((__int128)(ak)) * delta[0])) + delta_minus1);
} else {
delta[k] = ((0 - (((__int128)(ak)) * delta[(k - 1)])) + delta[(k - 2)]);
}
len = (k + 1);
if (qk > B) {
extra = (extra - 1);
if (extra <= 0) {
break;
}
}
k = (k + 1);
if (k > 500) {
break;
}
}
int32_t max_i = (len - 1);
int32_t* b_digits = (int32_t*)(((int32_t*)(calloc(((int64_t)((max_i + 1))), 4))));
__int128 beta_rem = beta_fp;
int32_t i = 1;
while (i <= max_i) {
__int128 d = delta[(i - 1)];
if (d == 0) {
b_digits[i] = 0;
} else {
__int128 bi = ceil_div_i128_i128_i128(beta_rem, d);
if (bi > ((__int128)(a[i]))) {
bi = ((__int128)(a[i]));
}
if (bi < 0) {
bi = 0;
}
b_digits[i] = ((int32_t)(bi));
beta_rem = ((bi * d) - beta_rem);
}
i = (i + 1);
}
int64_t* prefix = (int64_t*)(((int64_t*)(calloc(((int64_t)((max_i + 1))), 8))));
int64_t s = 0;
i = 1;
while (i <= max_i) {
s = (s + (((int64_t)(b_digits[i])) * q[(i - 1)]));
prefix[i] = s;
i = (i + 1);
}
int64_t* cand_r = (int64_t*)(((int64_t*)(calloc(200000, 8))));
int64_t* cand_l = (int64_t*)(((int64_t*)(calloc(200000, 8))));
int32_t nr = 0;
int32_t nl = 0;
cand_r[nr] = 0;
nr = (nr + 1);
cand_l[nl] = 0;
nl = (nl + 1);
int32_t k2 = 1;
while ((2 * k2) < (max_i + 1)) {
int32_t ie = (2 * k2);
int32_t io = ((2 * k2) - 1);
if (ie > max_i) {
break;
}
int64_t P = prefix[io];
int64_t st = q[io];
int32_t j = 0;
while (j < b_digits[ie]) {
int64_t n = (P + (((int64_t)(j)) * st));
if ((0 <= n && n <= B)) {
cand_r[nr] = n;
nr = (nr + 1);
}
j = (j + 1);
}
k2 = (k2 + 1);
}
k2 = 0;
while (((2 * k2) + 1) <= max_i) {
int32_t idx = (2 * k2);
int32_t inn = (idx + 1);
if (inn > max_i) {
break;
}
int64_t P = prefix[idx];
int64_t st = q[idx];
int32_t j = 0;
while (j < b_digits[inn]) {
int64_t n = (P + (((int64_t)(j)) * st));
if ((0 <= n && n <= B)) {
cand_l[nl] = n;
nl = (nl + 1);
}
j = (j + 1);
}
k2 = (k2 + 1);
}
__int128 best_r_gap = one();
__int128 best_l_gap = one();
int64_t best_r_n = 0;
int64_t best_l_n = 0;
i = 0;
while (i < nr) {
int64_t n = cand_r[i];
__int128 x = frac_mul_i128_i64(alpha_fp, n);
__int128 gap = (x - beta_fp);
if (gap < 0) {
gap = (gap + one());
}
if (gap < best_r_gap) {
best_r_gap = gap;
best_r_n = n;
}
i = (i + 1);
}
i = 0;
while (i < nl) {
int64_t n = cand_l[i];
__int128 x = frac_mul_i128_i64(alpha_fp, n);
__int128 gap = (beta_fp - x);
if (gap < 0) {
gap = (gap + one());
}
if (gap < best_l_gap) {
best_l_gap = gap;
best_l_n = n;
}
i = (i + 1);
}
int64_t ans = 0;
if (best_r_gap < best_l_gap) {
ans = best_r_n;
} else {
ans = best_l_n;
}
free(((void*)(delta)));
free(((void*)(a)));
free(((void*)(q)));
free(((void*)(b_digits)));
free(((void*)(prefix)));
free(((void*)(cand_r)));
free(((void*)(cand_l)));
return ans;
}
int64_t bqa_pi_d_i32_i64_i128_i128(int32_t d, int64_t n, __int128 pi_fp, __int128 beta_fp) {
int32_t* a0_ptr = (int32_t*)(((int32_t*)(calloc(1, 4))));
int32_t* period = (int32_t*)(((int32_t*)(calloc(4096, 4))));
int32_t* plen_ptr = (int32_t*)(((int32_t*)(calloc(1, 4))));
sqrt_cf_period_i32_ptr_i32_ptr_i32_ptr_i32(d, a0_ptr, period, plen_ptr);
int32_t a0 = a0_ptr[0];
int32_t plen = plen_ptr[0];
__int128 sqrt_D_fp = isqrt_scaled_i32(d);
__int128 alpha_fp = (sqrt_D_fp - FLOW_CHECKED_SHL((((__int128)(a0))), (96)));
double sqrt_d = sqrt(((double)(d)));
double pi_d = 3.14159265358979323846;
int64_t Bpos = ((int64_t)(floor(((pi_d + ((double)(n))) / sqrt_d))));
if (Bpos < 0) {
Bpos = 0;
}
if (Bpos > n) {
Bpos = n;
}
int64_t Bneg = ((int64_t)(floor(((((double)(n)) - pi_d) / sqrt_d))));
if (Bneg < 0) {
Bneg = 0;
}
if (Bneg > n) {
Bneg = n;
}
int64_t b_pos = best_b_positive_i128_i128_i64_ptr_i32_i32(alpha_fp, beta_fp, Bpos, period, plen);
__int128 one_m_beta = (one() - beta_fp);
int64_t t = best_b_positive_i128_i128_i64_ptr_i32_i32(alpha_fp, one_m_beta, Bneg, period, plen);
int64_t b_neg = (0 - t);
__int128 err1_fp = compute_error_fp_i128_i128_i64(sqrt_D_fp, pi_fp, b_pos);
__int128 err2_fp = compute_error_fp_i128_i128_i64(sqrt_D_fp, pi_fp, b_neg);
int64_t b_chosen = b_pos;
if (err2_fp < err1_fp) {
b_chosen = b_neg;
}
double a_real = (pi_d - (sqrt_d * ((double)(b_chosen))));
int64_t a_out = ((int64_t)(round(a_real)));
if (a_out > n) {
a_out = n;
}
if (a_out < (0 - n)) {
a_out = (0 - n);
}
free(((void*)(a0_ptr)));
free(((void*)(period)));
free(((void*)(plen_ptr)));
return a_out;
}
int32_t main(void) {
__int128 pi_fp = compute_pi_fp();
__int128 beta_fp = (pi_fp - FLOW_CHECKED_SHL((((__int128)(3))), (96)));
int64_t total = 0;
int64_t n = 10000000000000;
int32_t d = 2;
while (d < 100) {
if ((!(is_square_i32(d)))) {
int64_t a = bqa_pi_d_i32_i64_i128_i128(d, n, pi_fp, beta_fp);
if (a >= 0) {
total = (total + a);
} else {
total = (total + (0 - a));
}
}
d = (d + 1);
}
printf("%lld\n", total);
return 0;
}