# Project Euler 864
# Square + 1 = Squarefree: count squarefree x^2+1 for 1<=x<=n.
# Ported from native C to pure Flow.
import euler.nt { isqrt }
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
const N_VAL: i64 = 123567101113
const D_VAL: i64 = 30000000
# Globals for roots_minus_one_mod_p2 output
let mut g_r1: i64 = 0
let mut g_r2: i64 = 0
# Globals for negative_pell output
let mut g_pell_x: i64 = 0
let mut g_pell_y: i64 = 0
# dms globals
let mut g_dms_n: i64 = 0
let mut g_dms_D: i64 = 0
let mut g_dms_ans: i64 = 0
let mut g_dms_primes: ptr<i32> = null
let mut g_dms_roots: ptr<i64> = null
let mut g_dms_P: i32 = 0
# Mobius sieve globals
let mut g_spf_kmax: ptr<i32> = null
let mut g_mu_kmax: ptr<i8> = null
# Factoring globals
let mut g_primes_for_fact: ptr<i32> = null
let mut g_primes_for_fact_count: i32 = 0
function mod_pow_u32(a0: i32, e0: i64, p: i32) -> i32 {
let mut r: i64 = 1
let mut b: i64 = (a0 as i64) % (p as i64)
let mut e: i64 = e0
while e > 0 {
if (e & 1) == 1 { r = r * b % (p as i64) }
b = b * b % (p as i64)
e = e >> 1
}
return r as i32
}
# Modular inverse of m mod p2 using extended GCD with i128
function mod_inv_i128(m: i64, p2: i64) -> i64 {
let mut old_r: i128 = (m as i128) % (p2 as i128)
let mut r: i128 = p2 as i128
let mut old_s: i128 = 1
let mut s: i128 = 0
while r != 0 {
let q: i128 = old_r / r
let tmp: i128 = old_r - q * r
old_r = r
r = tmp
let tmp2: i128 = old_s - q * s
old_s = s
s = tmp2
}
return ((old_s % (p2 as i128) + (p2 as i128)) % (p2 as i128)) as i64
}
function primes_1mod4_upto(limit: i64) -> ptr<i32> {
if limit < 5 { return null }
let size: i32 = (limit / 2 + 1) as i32
let sieve: ptr<i8> = calloc(size as i64, 1)
let mut i: i32 = 0
while i < size {
sieve[i] = 1
i = i + 1
}
sieve[0] = 0
let r: i32 = isqrt(limit) as i32
let mut p: i32 = 3
while p <= r {
if sieve[p / 2] != 0 {
let mut j: i32 = (p * p) / 2
while j < size {
sieve[j] = 0
j = j + p
}
}
p = p + 2
}
let mut cnt: i32 = 0
i = 1
while i < size {
if sieve[i] != 0 && (2 * i + 1) % 4 == 1 { cnt = cnt + 1 }
i = i + 1
}
let result: ptr<i32> = calloc((cnt + 1) as i64, 4)
let mut idx: i32 = 0
i = 1
while i < size {
if sieve[i] != 0 && (2 * i + 1) % 4 == 1 {
result[idx] = 2 * i + 1
idx = idx + 1
}
i = i + 1
}
free(sieve)
return result
}
function sqrt_minus_one_mod_p(p: i32) -> i32 {
let exp: i64 = ((p - 1) as i64) / 2
let mut g: i32 = 2
while mod_pow_u32(g, exp, p) != p - 1 {
g = g + 1
}
return mod_pow_u32(g, ((p - 1) as i64) / 4, p)
}
function roots_minus_one_mod_p2(p: i32) -> void {
let r: i32 = sqrt_minus_one_mod_p(p)
let p2: i64 = (p as i64) * (p as i64)
let s: i64 = ((r as i64) * (r as i64) + 1) / (p as i64)
let inv: i32 = mod_pow_u32((2 * r) % p, (p - 2) as i64, p)
let t: i64 = (((p - (s % (p as i64))) as i64) * (inv as i64)) % (p as i64)
let R: i64 = ((r as i64) + t * (p as i64)) % p2
g_r1 = R
g_r2 = (p2 - R) % p2
}
function crt_combine(out: ptr<i64>, residues: ptr<i64>, len: i32, m: i64, a1: i64, a2: i64, p2: i64) -> void {
let inv: i64 = mod_inv_i128(m, p2)
let mut i: i32 = 0
while i < len {
let r: i64 = residues[i]
let r_mod: i64 = r % p2
let diff1: i64 = (a1 + p2 - r_mod) % p2
let t1: i64 = ((diff1 as i128) * (inv as i128) % (p2 as i128)) as i64
out[2 * i] = (r + m * t1) as i64
let diff2: i64 = (a2 + p2 - r_mod) % p2
let t2: i64 = ((diff2 as i128) * (inv as i128) % (p2 as i128)) as i64
out[2 * i + 1] = (r + m * t2) as i64
i = i + 1
}
}
function dms_dfs(start_idx: i32, d: i64, mod: i64, residues: ptr<i64>, len: i32, mu_sign: i32) -> void {
let mut i: i32 = start_idx
while i < g_dms_P {
let p: i32 = g_dms_primes[i]
let nd: i64 = d * (p as i64)
if nd > g_dms_D { break }
let p2: i64 = (p as i64) * (p as i64)
let a1: i64 = g_dms_roots[2 * i]
let a2: i64 = g_dms_roots[2 * i + 1]
let nmod: i64 = ((mod as i128) * (p2 as i128)) as i64
let nlen: i32 = len * 2
let nres: ptr<i64> = calloc(nlen as i64, 8)
crt_combine(nres, residues, len, mod, a1, a2, p2)
let mut A: i64 = 0
if nmod <= g_dms_n {
let q: i64 = g_dms_n / nmod
let rem: i64 = g_dms_n % nmod
A = q * (nlen as i64)
let mut j: i32 = 0
while j < nlen {
if nres[j] <= rem { A = A + 1 }
j = j + 1
}
} else {
let mut j: i32 = 0
while j < nlen {
if nres[j] <= g_dms_n { A = A + 1 }
j = j + 1
}
}
let neg_mu: i64 = (0 - (mu_sign as i64)) * A
g_dms_ans = g_dms_ans + neg_mu
dms_dfs(i + 1, nd, nmod, nres, nlen, 0 - mu_sign)
free(nres)
i = i + 1
}
}
function direct_mobius_sum(n: i64, D: i64) -> i64 {
let primes: ptr<i32> = primes_1mod4_upto(D)
# Count primes
let mut pcnt: i32 = 0
let mut pi: i32 = 0
while true {
let p: i32 = primes[pi]
if p == 0 { break }
pcnt = pcnt + 1
pi = pi + 1
}
g_dms_roots = calloc((pcnt * 2) as i64, 8)
let mut i: i32 = 0
while i < pcnt {
roots_minus_one_mod_p2(primes[i])
g_dms_roots[2 * i] = g_r1
g_dms_roots[2 * i + 1] = g_r2
i = i + 1
}
g_dms_n = n
g_dms_D = D
g_dms_primes = primes
g_dms_P = pcnt
g_dms_ans = n
let root_res: ptr<i64> = calloc(1, 8)
root_res[0] = 0
dms_dfs(0, 1, 1, root_res, 1, 1)
free(root_res)
free(g_dms_roots)
free(primes)
return g_dms_ans
}
function negative_pell_fundamental(D: i64, x_limit: i64) -> i32 {
let a0: i64 = isqrt(D)
if a0 * a0 == D { return 0 }
let mut m: i64 = 0
let mut d: i64 = 1
let mut a: i64 = a0
let mut p_prev: i64 = 1
let mut p: i64 = a0
let mut q_prev: i64 = 0
let mut q: i64 = 1
let mut period: i32 = 0
while true {
m = d * a - m
d = (D - m * m) / d
a = (a0 + m) / d
let new_p: i64 = a * p + p_prev
let new_q: i64 = a * q + q_prev
p_prev = p
p = new_p
q_prev = q
q = new_q
period = period + 1
if p_prev > x_limit { return 0 }
if a == 2 * a0 { break }
}
if period % 2 == 0 { return 0 }
if p_prev * p_prev - D * q_prev * q_prev != 0 - 1 { return 0 }
g_pell_x = p_prev
g_pell_y = q_prev
return 1
}
function linear_sieve_kmax(Kmax: i64) -> void {
g_spf_kmax = calloc(Kmax + 1, 4)
g_mu_kmax = calloc(Kmax + 1, 1)
let primes: ptr<i32> = calloc(Kmax / 10 + 10000, 4)
let mut pc: i32 = 0
g_mu_kmax[1] = 1
let mut i: i32 = 2
while (i as i64) <= Kmax {
if g_spf_kmax[i] == 0 {
g_spf_kmax[i] = i
primes[pc] = i
pc = pc + 1
g_mu_kmax[i] = 0 - 1
}
let mut j: i32 = 0
while j < pc {
let p: i32 = primes[j]
let ip: i64 = (i as i64) * (p as i64)
if ip > Kmax { break }
g_spf_kmax[ip as i32] = p
if i % p == 0 {
g_mu_kmax[ip as i32] = 0
break
}
g_mu_kmax[ip as i32] = 0 - g_mu_kmax[i]
j = j + 1
}
i = i + 1
}
free(primes)
}
function build_primes_for_factoring(limit: i64) -> void {
g_primes_for_fact = null
g_primes_for_fact_count = 0
if limit < 2 { return }
let size: i32 = (limit / 2 + 1) as i32
let sieve: ptr<i8> = calloc(size as i64, 1)
let mut i: i32 = 0
while i < size {
sieve[i] = 1
i = i + 1
}
sieve[0] = 0
let r: i32 = isqrt(limit) as i32
let mut p: i32 = 3
while p <= r {
if sieve[p / 2] != 0 {
let mut j: i32 = (p * p) / 2
while j < size {
sieve[j] = 0
j = j + p
}
}
p = p + 2
}
let mut cnt: i32 = 1
i = 1
while i < size {
if sieve[i] != 0 { cnt = cnt + 1 }
i = i + 1
}
g_primes_for_fact = calloc(cnt as i64, 4)
g_primes_for_fact[g_primes_for_fact_count] = 2
g_primes_for_fact_count = g_primes_for_fact_count + 1
i = 1
while i < size {
if sieve[i] != 0 {
g_primes_for_fact[g_primes_for_fact_count] = 2 * i + 1
g_primes_for_fact_count = g_primes_for_fact_count + 1
}
i = i + 1
}
free(sieve)
}
function factor_distinct_primes(n: i64, out: ptr<i32>) -> i32 {
let mut cnt: i32 = 0
let mut t: i64 = n
let mut i: i32 = 0
while i < g_primes_for_fact_count {
let p: i32 = g_primes_for_fact[i]
if (p as i64) * (p as i64) > t { break }
if t % (p as i64) == 0 {
out[cnt] = p
cnt = cnt + 1
while t % (p as i64) == 0 { t = t / (p as i64) }
}
i = i + 1
}
if t > 1 {
out[cnt] = t as i32
cnt = cnt + 1
}
return cnt
}
function mobius_tail_sum_for_y(y: i64, D: i64) -> i64 {
let pf: ptr<i32> = calloc(64, 4)
let npf: i32 = factor_distinct_primes(y, pf)
let mut total: i64 = 0
let subsets: i32 = 1 << npf
let mut mask: i32 = 0
while mask < subsets {
let mut prod: i64 = 1
let mut parity: i32 = 0
let mut b: i32 = 0
while b < npf {
if (mask & (1 << b)) != 0 {
prod = prod * (pf[b] as i64)
parity = parity ^ 1
}
b = b + 1
}
if prod > D {
if parity != 0 {
total = total - 1
} else {
total = total + 1
}
}
mask = mask + 1
}
free(pf)
return total
}
function correction_via_pell(n: i64, D: i64) -> i64 {
let Kmax: i64 = (((n as i128) * (n as i128) + 1) / ((D as i128) * (D as i128))) as i64
if Kmax < 2 { return 0 }
linear_sieve_kmax(Kmax)
let isqrt_n: i64 = isqrt(n)
build_primes_for_factoring(isqrt_n + 1)
let mut corr: i64 = 0
let mut k: i64 = 2
while k <= Kmax {
if g_mu_kmax[k] != 0 {
# Filter: if an odd prime p == 3 (mod 4) divides k, skip
let mut t: i64 = k
let mut ok: i32 = 1
while t > 1 {
let p: i32 = g_spf_kmax[t as i32]
t = t / (p as i64)
if p != 2 && (p & 3) == 3 {
ok = 0
break
}
}
if ok != 0 {
if negative_pell_fundamental(k, n) != 0 {
let x: i64 = g_pell_x
let y: i64 = g_pell_y
let A: i128 = (x as i128) * (x as i128) + (k as i128) * (y as i128) * (y as i128)
let B: i128 = (2 as i128) * (x as i128) * (y as i128)
let mut cx: i64 = x
let mut cy: i64 = y
while cx <= n {
if cy > D {
corr = corr + mobius_tail_sum_for_y(cy, D)
}
let new_x: i128 = A * (cx as i128) + (k as i128) * B * (cy as i128)
let new_y: i128 = B * (cx as i128) + A * (cy as i128)
if new_x > (n as i128) { break }
cx = new_x as i64
cy = new_y as i64
}
}
}
}
k = k + 1
}
free(g_spf_kmax)
free(g_mu_kmax)
free(g_primes_for_fact)
return corr
}
function main() -> i32 {
let direct: i64 = direct_mobius_sum(N_VAL, D_VAL)
let corr: i64 = correction_via_pell(N_VAL, D_VAL)
printf("%lld\n", direct + corr)
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);
int32_t mod_pow_u32_i32_i64_i32(int32_t a0, int64_t e0, int32_t p);
int64_t mod_inv_i128_i64_i64(int64_t m, int64_t p2);
int32_t* primes_1mod4_upto_i64(int64_t limit);
int32_t sqrt_minus_one_mod_p_i32(int32_t p);
void roots_minus_one_mod_p2_i32(int32_t p);
void crt_combine_ptr_i64_ptr_i64_i32_i64_i64_i64_i64(int64_t* out, int64_t* residues, int32_t len, int64_t m, int64_t a1, int64_t a2, int64_t p2);
void dms_dfs_i32_i64_i64_ptr_i64_i32_i32(int32_t start_idx, int64_t d, int64_t mod, int64_t* residues, int32_t len, int32_t mu_sign);
int64_t direct_mobius_sum_i64_i64(int64_t n, int64_t D);
int32_t negative_pell_fundamental_i64_i64(int64_t D, int64_t x_limit);
void linear_sieve_kmax_i64(int64_t Kmax);
void build_primes_for_factoring_i64(int64_t limit);
int32_t factor_distinct_primes_i64_ptr_i32(int64_t n, int32_t* out);
int64_t mobius_tail_sum_for_y_i64_i64(int64_t y, int64_t D);
int64_t correction_via_pell_i64_i64(int64_t n, int64_t D);
int32_t main(void);
static const int64_t N_VAL = 123567101113;
static const int64_t D_VAL = 30000000;
/* Module statics */
static int64_t g_r1 = 0;
static int64_t g_r2 = 0;
static int64_t g_pell_x = 0;
static int64_t g_pell_y = 0;
static int64_t g_dms_n = 0;
static int64_t g_dms_D = 0;
static int64_t g_dms_ans = 0;
static int32_t* g_dms_primes = NULL;
static int64_t* g_dms_roots = NULL;
static int32_t g_dms_P = 0;
static int32_t* g_spf_kmax = NULL;
static int8_t* g_mu_kmax = NULL;
static int32_t* g_primes_for_fact = NULL;
static int32_t g_primes_for_fact_count = 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;
}
int32_t mod_pow_u32_i32_i64_i32(int32_t a0, int64_t e0, int32_t p) {
int64_t r = 1;
int64_t b = FLOW_CHECKED_MOD((((int64_t)(a0))), (((int64_t)(p))));
int64_t e = e0;
while (e > 0) {
if ((e & 1) == 1) {
r = FLOW_CHECKED_MOD(((r * b)), (((int64_t)(p))));
}
b = FLOW_CHECKED_MOD(((b * b)), (((int64_t)(p))));
e = FLOW_CHECKED_SHR((e), (1));
}
return ((int32_t)(r));
}
int64_t mod_inv_i128_i64_i64(int64_t m, int64_t p2) {
__int128 old_r = FLOW_CHECKED_MOD((((__int128)(m))), (((__int128)(p2))));
__int128 r = ((__int128)(p2));
__int128 old_s = 1;
__int128 s = 0;
while (r != 0) {
__int128 q = FLOW_CHECKED_DIV((old_r), (r));
__int128 tmp = (old_r - (q * r));
old_r = r;
r = tmp;
__int128 tmp2 = (old_s - (q * s));
old_s = s;
s = tmp2;
}
return ((int64_t)(FLOW_CHECKED_MOD(((FLOW_CHECKED_MOD((old_s), (((__int128)(p2)))) + ((__int128)(p2)))), (((__int128)(p2))))));
}
int32_t* primes_1mod4_upto_i64(int64_t limit) {
if (limit < 5) {
return NULL;
}
int32_t size = ((int32_t)((FLOW_CHECKED_DIV((limit), (2)) + 1)));
int8_t* sieve = (int8_t*)(calloc(((int64_t)(size)), 1));
int32_t i = 0;
while (i < size) {
sieve[i] = 1;
i = (i + 1);
}
sieve[0] = 0;
int32_t r = ((int32_t)(isqrt_i64(limit)));
int32_t p = 3;
while (p <= r) {
if (sieve[FLOW_CHECKED_DIV((p), (2))] != 0) {
int32_t j = FLOW_CHECKED_DIV(((p * p)), (2));
while (j < size) {
sieve[j] = 0;
j = (j + p);
}
}
p = (p + 2);
}
int32_t cnt = 0;
i = 1;
while (i < size) {
if ((sieve[i] != 0 && FLOW_CHECKED_MOD((((2 * i) + 1)), (4)) == 1)) {
cnt = (cnt + 1);
}
i = (i + 1);
}
int32_t* result = (int32_t*)(calloc(((int64_t)((cnt + 1))), 4));
int32_t idx = 0;
i = 1;
while (i < size) {
if ((sieve[i] != 0 && FLOW_CHECKED_MOD((((2 * i) + 1)), (4)) == 1)) {
result[idx] = ((2 * i) + 1);
idx = (idx + 1);
}
i = (i + 1);
}
free(sieve);
return result;
}
int32_t sqrt_minus_one_mod_p_i32(int32_t p) {
int64_t exp = FLOW_CHECKED_DIV((((int64_t)((p - 1)))), (2));
int32_t g = 2;
while (mod_pow_u32_i32_i64_i32(g, exp, p) != (p - 1)) {
g = (g + 1);
}
return mod_pow_u32_i32_i64_i32(g, FLOW_CHECKED_DIV((((int64_t)((p - 1)))), (4)), p);
}
void roots_minus_one_mod_p2_i32(int32_t p) {
int32_t r = sqrt_minus_one_mod_p_i32(p);
int64_t p2 = (((int64_t)(p)) * ((int64_t)(p)));
int64_t s = FLOW_CHECKED_DIV((((((int64_t)(r)) * ((int64_t)(r))) + 1)), (((int64_t)(p))));
int32_t inv = mod_pow_u32_i32_i64_i32(FLOW_CHECKED_MOD(((2 * r)), (p)), ((int64_t)((p - 2))), p);
int64_t t = FLOW_CHECKED_MOD(((((int64_t)((p - FLOW_CHECKED_MOD((s), (((int64_t)(p))))))) * ((int64_t)(inv)))), (((int64_t)(p))));
int64_t R = FLOW_CHECKED_MOD(((((int64_t)(r)) + (t * ((int64_t)(p))))), (p2));
g_r1 = R;
g_r2 = FLOW_CHECKED_MOD(((p2 - R)), (p2));
}
void crt_combine_ptr_i64_ptr_i64_i32_i64_i64_i64_i64(int64_t* out, int64_t* residues, int32_t len, int64_t m, int64_t a1, int64_t a2, int64_t p2) {
int64_t inv = mod_inv_i128_i64_i64(m, p2);
int32_t i = 0;
while (i < len) {
int64_t r = residues[i];
int64_t r_mod = FLOW_CHECKED_MOD((r), (p2));
int64_t diff1 = FLOW_CHECKED_MOD((((a1 + p2) - r_mod)), (p2));
int64_t t1 = ((int64_t)(FLOW_CHECKED_MOD(((((__int128)(diff1)) * ((__int128)(inv)))), (((__int128)(p2))))));
out[(2 * i)] = ((int64_t)((r + (m * t1))));
int64_t diff2 = FLOW_CHECKED_MOD((((a2 + p2) - r_mod)), (p2));
int64_t t2 = ((int64_t)(FLOW_CHECKED_MOD(((((__int128)(diff2)) * ((__int128)(inv)))), (((__int128)(p2))))));
out[((2 * i) + 1)] = ((int64_t)((r + (m * t2))));
i = (i + 1);
}
}
void dms_dfs_i32_i64_i64_ptr_i64_i32_i32(int32_t start_idx, int64_t d, int64_t mod, int64_t* residues, int32_t len, int32_t mu_sign) {
int32_t i = start_idx;
while (i < g_dms_P) {
int32_t p = g_dms_primes[i];
int64_t nd = (d * ((int64_t)(p)));
if (nd > g_dms_D) {
break;
}
int64_t p2 = (((int64_t)(p)) * ((int64_t)(p)));
int64_t a1 = g_dms_roots[(2 * i)];
int64_t a2 = g_dms_roots[((2 * i) + 1)];
int64_t nmod = ((int64_t)((((__int128)(mod)) * ((__int128)(p2)))));
int32_t nlen = (len * 2);
int64_t* nres = (int64_t*)(calloc(((int64_t)(nlen)), 8));
crt_combine_ptr_i64_ptr_i64_i32_i64_i64_i64_i64(nres, residues, len, mod, a1, a2, p2);
int64_t A = 0;
if (nmod <= g_dms_n) {
int64_t q = FLOW_CHECKED_DIV((g_dms_n), (nmod));
int64_t rem = FLOW_CHECKED_MOD((g_dms_n), (nmod));
A = (q * ((int64_t)(nlen)));
int32_t j = 0;
while (j < nlen) {
if (nres[j] <= rem) {
A = (A + 1);
}
j = (j + 1);
}
} else {
int32_t j = 0;
while (j < nlen) {
if (nres[j] <= g_dms_n) {
A = (A + 1);
}
j = (j + 1);
}
}
int64_t neg_mu = ((0 - ((int64_t)(mu_sign))) * A);
g_dms_ans = (g_dms_ans + neg_mu);
dms_dfs_i32_i64_i64_ptr_i64_i32_i32((i + 1), nd, nmod, nres, nlen, (0 - mu_sign));
free(nres);
i = (i + 1);
}
}
int64_t direct_mobius_sum_i64_i64(int64_t n, int64_t D) {
int32_t* primes = (int32_t*)(primes_1mod4_upto_i64(D));
int32_t pcnt = 0;
int32_t pi = 0;
while (1) {
int32_t p = primes[pi];
if (p == 0) {
break;
}
pcnt = (pcnt + 1);
pi = (pi + 1);
}
g_dms_roots = calloc(((int64_t)((pcnt * 2))), 8);
int32_t i = 0;
while (i < pcnt) {
roots_minus_one_mod_p2_i32(primes[i]);
g_dms_roots[(2 * i)] = g_r1;
g_dms_roots[((2 * i) + 1)] = g_r2;
i = (i + 1);
}
g_dms_n = n;
g_dms_D = D;
g_dms_primes = primes;
g_dms_P = pcnt;
g_dms_ans = n;
int64_t* root_res = (int64_t*)(calloc(1, 8));
root_res[0] = 0;
dms_dfs_i32_i64_i64_ptr_i64_i32_i32(0, 1, 1, root_res, 1, 1);
free(root_res);
free(g_dms_roots);
free(primes);
return g_dms_ans;
}
int32_t negative_pell_fundamental_i64_i64(int64_t D, int64_t x_limit) {
int64_t a0 = isqrt_i64(D);
if ((a0 * a0) == D) {
return 0;
}
int64_t m = 0;
int64_t d = 1;
int64_t a = a0;
int64_t p_prev = 1;
int64_t p = a0;
int64_t q_prev = 0;
int64_t q = 1;
int32_t period = 0;
while (1) {
m = ((d * a) - m);
d = FLOW_CHECKED_DIV(((D - (m * m))), (d));
a = FLOW_CHECKED_DIV(((a0 + m)), (d));
int64_t new_p = ((a * p) + p_prev);
int64_t new_q = ((a * q) + q_prev);
p_prev = p;
p = new_p;
q_prev = q;
q = new_q;
period = (period + 1);
if (p_prev > x_limit) {
return 0;
}
if (a == (2 * a0)) {
break;
}
}
if (FLOW_CHECKED_MOD((period), (2)) == 0) {
return 0;
}
if (((p_prev * p_prev) - ((D * q_prev) * q_prev)) != (0 - 1)) {
return 0;
}
g_pell_x = p_prev;
g_pell_y = q_prev;
return 1;
}
void linear_sieve_kmax_i64(int64_t Kmax) {
g_spf_kmax = calloc((Kmax + 1), 4);
g_mu_kmax = calloc((Kmax + 1), 1);
int32_t* primes = (int32_t*)(calloc((FLOW_CHECKED_DIV((Kmax), (10)) + 10000), 4));
int32_t pc = 0;
g_mu_kmax[1] = 1;
int32_t i = 2;
while (((int64_t)(i)) <= Kmax) {
if (g_spf_kmax[i] == 0) {
g_spf_kmax[i] = i;
primes[pc] = i;
pc = (pc + 1);
g_mu_kmax[i] = (0 - 1);
}
int32_t j = 0;
while (j < pc) {
int32_t p = primes[j];
int64_t ip = (((int64_t)(i)) * ((int64_t)(p)));
if (ip > Kmax) {
break;
}
g_spf_kmax[((int32_t)(ip))] = p;
if (FLOW_CHECKED_MOD((i), (p)) == 0) {
g_mu_kmax[((int32_t)(ip))] = 0;
break;
}
g_mu_kmax[((int32_t)(ip))] = (0 - g_mu_kmax[i]);
j = (j + 1);
}
i = (i + 1);
}
free(primes);
}
void build_primes_for_factoring_i64(int64_t limit) {
g_primes_for_fact = NULL;
g_primes_for_fact_count = 0;
if (limit < 2) {
return;
}
int32_t size = ((int32_t)((FLOW_CHECKED_DIV((limit), (2)) + 1)));
int8_t* sieve = (int8_t*)(calloc(((int64_t)(size)), 1));
int32_t i = 0;
while (i < size) {
sieve[i] = 1;
i = (i + 1);
}
sieve[0] = 0;
int32_t r = ((int32_t)(isqrt_i64(limit)));
int32_t p = 3;
while (p <= r) {
if (sieve[FLOW_CHECKED_DIV((p), (2))] != 0) {
int32_t j = FLOW_CHECKED_DIV(((p * p)), (2));
while (j < size) {
sieve[j] = 0;
j = (j + p);
}
}
p = (p + 2);
}
int32_t cnt = 1;
i = 1;
while (i < size) {
if (sieve[i] != 0) {
cnt = (cnt + 1);
}
i = (i + 1);
}
g_primes_for_fact = calloc(((int64_t)(cnt)), 4);
g_primes_for_fact[g_primes_for_fact_count] = 2;
g_primes_for_fact_count = (g_primes_for_fact_count + 1);
i = 1;
while (i < size) {
if (sieve[i] != 0) {
g_primes_for_fact[g_primes_for_fact_count] = ((2 * i) + 1);
g_primes_for_fact_count = (g_primes_for_fact_count + 1);
}
i = (i + 1);
}
free(sieve);
}
int32_t factor_distinct_primes_i64_ptr_i32(int64_t n, int32_t* out) {
int32_t cnt = 0;
int64_t t = n;
int32_t i = 0;
while (i < g_primes_for_fact_count) {
int32_t p = g_primes_for_fact[i];
if ((((int64_t)(p)) * ((int64_t)(p))) > t) {
break;
}
if (FLOW_CHECKED_MOD((t), (((int64_t)(p)))) == 0) {
out[cnt] = p;
cnt = (cnt + 1);
while (FLOW_CHECKED_MOD((t), (((int64_t)(p)))) == 0) {
t = FLOW_CHECKED_DIV((t), (((int64_t)(p))));
}
}
i = (i + 1);
}
if (t > 1) {
out[cnt] = ((int32_t)(t));
cnt = (cnt + 1);
}
return cnt;
}
int64_t mobius_tail_sum_for_y_i64_i64(int64_t y, int64_t D) {
int32_t* pf = (int32_t*)(calloc(64, 4));
int32_t npf = factor_distinct_primes_i64_ptr_i32(y, pf);
int64_t total = 0;
int32_t subsets = FLOW_CHECKED_SHL((1), (npf));
int32_t mask = 0;
while (mask < subsets) {
int64_t prod = 1;
int32_t parity = 0;
int32_t b = 0;
while (b < npf) {
if ((mask & FLOW_CHECKED_SHL((1), (b))) != 0) {
prod = (prod * ((int64_t)(pf[b])));
parity = (parity ^ 1);
}
b = (b + 1);
}
if (prod > D) {
if (parity != 0) {
total = (total - 1);
} else {
total = (total + 1);
}
}
mask = (mask + 1);
}
free(pf);
return total;
}
int64_t correction_via_pell_i64_i64(int64_t n, int64_t D) {
int64_t Kmax = ((int64_t)(FLOW_CHECKED_DIV((((((__int128)(n)) * ((__int128)(n))) + 1)), ((((__int128)(D)) * ((__int128)(D)))))));
if (Kmax < 2) {
return 0;
}
linear_sieve_kmax_i64(Kmax);
int64_t isqrt_n = isqrt_i64(n);
build_primes_for_factoring_i64((isqrt_n + 1));
int64_t corr = 0;
int64_t k = 2;
while (k <= Kmax) {
if (g_mu_kmax[k] != 0) {
int64_t t = k;
int32_t ok = 1;
while (t > 1) {
int32_t p = g_spf_kmax[((int32_t)(t))];
t = FLOW_CHECKED_DIV((t), (((int64_t)(p))));
if ((p != 2 && (p & 3) == 3)) {
ok = 0;
break;
}
}
if (ok != 0) {
if (negative_pell_fundamental_i64_i64(k, n) != 0) {
int64_t x = g_pell_x;
int64_t y = g_pell_y;
__int128 A = ((((__int128)(x)) * ((__int128)(x))) + ((((__int128)(k)) * ((__int128)(y))) * ((__int128)(y))));
__int128 B = ((((__int128)(2)) * ((__int128)(x))) * ((__int128)(y)));
int64_t cx = x;
int64_t cy = y;
while (cx <= n) {
if (cy > D) {
corr = (corr + mobius_tail_sum_for_y_i64_i64(cy, D));
}
__int128 new_x = ((A * ((__int128)(cx))) + ((((__int128)(k)) * B) * ((__int128)(cy))));
__int128 new_y = ((B * ((__int128)(cx))) + (A * ((__int128)(cy))));
if (new_x > ((__int128)(n))) {
break;
}
cx = ((int64_t)(new_x));
cy = ((int64_t)(new_y));
}
}
}
}
k = (k + 1);
}
free(g_spf_kmax);
free(g_mu_kmax);
free(g_primes_for_fact);
return corr;
}
int32_t main(void) {
int64_t direct = direct_mobius_sum_i64_i64(N_VAL, D_VAL);
int64_t corr = correction_via_pell_i64_i64(N_VAL, D_VAL);
printf("%lld\n", (direct + corr));
return 0;
}