# Project Euler 372
# R(2e6, 1e9): count pairs with floor((y/x)^2) odd.
import euler.nt { isqrt }
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
}
function igcd(a0: i64, b0: i64) -> i64 {
let mut a: i64 = a0
let mut b: i64 = b0
if a < 0 { a = 0 - a }
if b < 0 { b = 0 - b }
while b != 0 {
let t: i64 = a % b
a = b
b = t
}
return a
}
function floor_sum(n0: i64, m0: i64, a0: i64, b0: i64) -> i64 {
let mut res: i64 = 0
let mut n: i64 = n0
let mut m: i64 = m0
let mut a: i64 = a0
let mut b: i64 = b0
while 1 == 1 {
if a >= m {
res = res + (n - 1) * n * (a / m) / 2
a = a % m
}
if b >= m {
res = res + n * (b / m)
b = b % m
}
let y_max: i64 = a * n + b
if y_max < m { break }
n = y_max / m
b = y_max % m
let tmp: i64 = a
a = m
m = tmp
}
return res
}
function isqrt128(n: i128) -> i64 {
if n <= (0 as i128) { return 0 }
if n < (2 as i128) { return n as i64 }
let mut x: i128 = n
let mut y: i128 = (x + (1 as i128)) / (2 as i128)
while y < x {
x = y
y = (x + n / x) / (2 as i128)
}
return x as i64
}
function floor_surd(a: i64, b: i64, c: i64, k: i64) -> i64 {
if b == 0 || k == 0 { return a / c }
let mut bb: i64 = b
if bb < 0 { bb = 0 - bb }
let bbk: i128 = (bb as i128) * (bb as i128) * (k as i128)
let s: i64 = isqrt128(bbk)
if b > 0 {
let t: i64 = (a + s) / c
let u: i64 = (t + 1) * c - a
if u <= 0 { return t + 1 }
if (u as i128) * (u as i128) <= bbk { return t + 1 }
return t
}
let t2: i64 = (a - s - 1) / c
let d: i64 = a - (t2 + 1) * c
if d >= 0 && (d as i128) * (d as i128) >= bbk { return t2 + 1 }
return t2
}
function normalize(a0: i64, b0: i64, c0: i64, out_a: ptr<i64>, out_b: ptr<i64>, out_c: ptr<i64>) -> void {
let mut a: i64 = a0
let mut b: i64 = b0
let mut c: i64 = c0
if c < 0 {
a = 0 - a
b = 0 - b
c = 0 - c
}
let g: i64 = igcd(igcd(a, b), c)
if g > 1 {
a = a / g
b = b / g
c = c / g
}
out_a[0] = a
out_b[0] = b
out_c[0] = c
}
function sum_floor_mul_surd(a0: i64, b0: i64, c0: i64, k: i64, n0: i64) -> i64 {
if n0 <= 0 { return 0 }
if b0 == 0 || k == 0 {
let p: i64 = a0
let q: i64 = c0
if p >= 0 {
return floor_sum(n0, q, p, p)
}
let pp: i64 = 0 - p
return 0 - floor_sum(n0, q, pp, pp + q - 1)
}
let mut aa: i64 = a0
let mut bb: i64 = b0
let mut cc: i64 = c0
let mut nn: i64 = n0
let mut res: i64 = 0
let mut sign: i64 = 1
let box_a: array<i64, 1> = [0]
let box_b: array<i64, 1> = [0]
let box_c: array<i64, 1> = [0]
while nn > 0 {
let q_int: i64 = floor_surd(aa, bb, cc, k)
if q_int != 0 {
res = res + sign * q_int * nn * (nn + 1) / 2
aa = aa - q_int * cc
}
if aa == 0 && bb == 0 { break }
let m: i64 = floor_surd(aa * nn, bb * nn, cc, k)
if m == 0 { break }
res = res + sign * nn * m
# reciprocal
let mut A: i64 = cc * aa
let mut B: i64 = 0 - cc * bb
let mut C: i64 = aa * aa - bb * bb * k
if C < 0 {
A = 0 - A
B = 0 - B
C = 0 - C
}
normalize(A, B, C, box_a, box_b, box_c)
aa = box_a[0]
bb = box_b[0]
cc = box_c[0]
nn = m
sign = 0 - sign
}
return res
}
function P_func(k: i64, L: i64, U: i64, len_side: i64, total_pairs: i64, U2m1: i64, Lm1: i64) -> i64 {
if k == 1 {
return len_side * (len_side - 1) / 2
}
let r: i64 = isqrt(k)
if r * r == k {
let s: i64 = r
let mut b: i64 = U / s
if b < L { return total_pairs }
if b > U { b = U }
let cnt: i64 = b - L + 1
let sum_x: i64 = (L + b) * cnt / 2
let partial: i64 = s * sum_x - L * cnt
let full_cnt: i64 = U - b
return partial + full_cnt * len_side
}
let mut b2: i64 = isqrt(U2m1 / k)
if b2 < L { return total_pairs }
if b2 > U { b2 = U }
let cnt2: i64 = b2 - L + 1
let sum_floor: i64 = sum_floor_mul_surd(0, 1, 1, k, b2) - sum_floor_mul_surd(0, 1, 1, k, Lm1)
let partial2: i64 = sum_floor - (L - 1) * cnt2
let full_cnt2: i64 = U - b2
return partial2 + full_cnt2 * len_side
}
function solve(M: i64, N: i64) -> i64 {
let L: i64 = M + 1
let U: i64 = N
if L > U { return 0 }
let len_side: i64 = U - L + 1
let total_pairs: i64 = len_side * len_side
let k_max: i64 = (U * U) / (L * L)
let mut last_odd: i64 = k_max
if (k_max & 1) == 0 { last_odd = k_max - 1 }
let U2m1: i64 = U * U - 1
let Lm1: i64 = L - 1
let mut ans: i64 = 0
let mut n: i64 = 1
while n <= last_odd {
ans = ans + P_func(n + 1, L, U, len_side, total_pairs, U2m1, Lm1) - P_func(n, L, U, len_side, total_pairs, U2m1, Lm1)
n = n + 2
}
return ans
}
function main() -> i32 {
let ans: i64 = solve(2000000, 1000000000)
printf("%lld\n", ans)
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);
int64_t igcd_i64_i64(int64_t a0, int64_t b0);
int64_t floor_sum_i64_i64_i64_i64(int64_t n0, int64_t m0, int64_t a0, int64_t b0);
int64_t isqrt128_i128(__int128 n);
int64_t floor_surd_i64_i64_i64_i64(int64_t a, int64_t b, int64_t c, int64_t k);
void normalize_i64_i64_i64_ptr_i64_ptr_i64_ptr_i64(int64_t a0, int64_t b0, int64_t c0, int64_t* out_a, int64_t* out_b, int64_t* out_c);
int64_t sum_floor_mul_surd_i64_i64_i64_i64_i64(int64_t a0, int64_t b0, int64_t c0, int64_t k, int64_t n0);
int64_t P_func_i64_i64_i64_i64_i64_i64_i64(int64_t k, int64_t L, int64_t U, int64_t len_side, int64_t total_pairs, int64_t U2m1, int64_t Lm1);
int64_t solve_i64_i64(int64_t M, int64_t N);
int32_t main(void);
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;
}
int64_t igcd_i64_i64(int64_t a0, int64_t b0) {
int64_t a = a0;
int64_t b = b0;
if (a < 0) {
a = (0 - a);
}
if (b < 0) {
b = (0 - b);
}
while (b != 0) {
int64_t t = FLOW_CHECKED_MOD((a), (b));
a = b;
b = t;
}
return a;
}
int64_t floor_sum_i64_i64_i64_i64(int64_t n0, int64_t m0, int64_t a0, int64_t b0) {
int64_t res = 0;
int64_t n = n0;
int64_t m = m0;
int64_t a = a0;
int64_t b = b0;
while (1 == 1) {
if (a >= m) {
res = (res + FLOW_CHECKED_DIV(((((n - 1) * n) * FLOW_CHECKED_DIV((a), (m)))), (2)));
a = FLOW_CHECKED_MOD((a), (m));
}
if (b >= m) {
res = (res + (n * FLOW_CHECKED_DIV((b), (m))));
b = FLOW_CHECKED_MOD((b), (m));
}
int64_t y_max = ((a * n) + b);
if (y_max < m) {
break;
}
n = FLOW_CHECKED_DIV((y_max), (m));
b = FLOW_CHECKED_MOD((y_max), (m));
int64_t tmp = a;
a = m;
m = tmp;
}
return res;
}
int64_t isqrt128_i128(__int128 n) {
if (n <= ((__int128)(0))) {
return 0;
}
if (n < ((__int128)(2))) {
return ((int64_t)(n));
}
__int128 x = n;
__int128 y = FLOW_CHECKED_DIV(((x + ((__int128)(1)))), (((__int128)(2))));
while (y < x) {
x = y;
y = FLOW_CHECKED_DIV(((x + FLOW_CHECKED_DIV((n), (x)))), (((__int128)(2))));
}
return ((int64_t)(x));
}
int64_t floor_surd_i64_i64_i64_i64(int64_t a, int64_t b, int64_t c, int64_t k) {
if ((b == 0 || k == 0)) {
return FLOW_CHECKED_DIV((a), (c));
}
int64_t bb = b;
if (bb < 0) {
bb = (0 - bb);
}
__int128 bbk = ((((__int128)(bb)) * ((__int128)(bb))) * ((__int128)(k)));
int64_t s = isqrt128_i128(bbk);
if (b > 0) {
int64_t t = FLOW_CHECKED_DIV(((a + s)), (c));
int64_t u = (((t + 1) * c) - a);
if (u <= 0) {
return (t + 1);
}
if ((((__int128)(u)) * ((__int128)(u))) <= bbk) {
return (t + 1);
}
return t;
}
int64_t t2 = FLOW_CHECKED_DIV((((a - s) - 1)), (c));
int64_t d = (a - ((t2 + 1) * c));
if ((d >= 0 && (((__int128)(d)) * ((__int128)(d))) >= bbk)) {
return (t2 + 1);
}
return t2;
}
void normalize_i64_i64_i64_ptr_i64_ptr_i64_ptr_i64(int64_t a0, int64_t b0, int64_t c0, int64_t* out_a, int64_t* out_b, int64_t* out_c) {
int64_t a = a0;
int64_t b = b0;
int64_t c = c0;
if (c < 0) {
a = (0 - a);
b = (0 - b);
c = (0 - c);
}
int64_t g = igcd_i64_i64(igcd_i64_i64(a, b), c);
if (g > 1) {
a = FLOW_CHECKED_DIV((a), (g));
b = FLOW_CHECKED_DIV((b), (g));
c = FLOW_CHECKED_DIV((c), (g));
}
out_a[0] = a;
out_b[0] = b;
out_c[0] = c;
}
int64_t sum_floor_mul_surd_i64_i64_i64_i64_i64(int64_t a0, int64_t b0, int64_t c0, int64_t k, int64_t n0) {
if (n0 <= 0) {
return 0;
}
if ((b0 == 0 || k == 0)) {
int64_t p = a0;
int64_t q = c0;
if (p >= 0) {
return floor_sum_i64_i64_i64_i64(n0, q, p, p);
}
int64_t pp = (0 - p);
return (0 - floor_sum_i64_i64_i64_i64(n0, q, pp, ((pp + q) - 1)));
}
int64_t aa = a0;
int64_t bb = b0;
int64_t cc = c0;
int64_t nn = n0;
int64_t res = 0;
int64_t sign = 1;
int64_t box_a[1] = { 0 };
int64_t box_b[1] = { 0 };
int64_t box_c[1] = { 0 };
while (nn > 0) {
int64_t q_int = floor_surd_i64_i64_i64_i64(aa, bb, cc, k);
if (q_int != 0) {
res = (res + FLOW_CHECKED_DIV(((((sign * q_int) * nn) * (nn + 1))), (2)));
aa = (aa - (q_int * cc));
}
if ((aa == 0 && bb == 0)) {
break;
}
int64_t m = floor_surd_i64_i64_i64_i64((aa * nn), (bb * nn), cc, k);
if (m == 0) {
break;
}
res = (res + ((sign * nn) * m));
int64_t A = (cc * aa);
int64_t B = (0 - (cc * bb));
int64_t C = ((aa * aa) - ((bb * bb) * k));
if (C < 0) {
A = (0 - A);
B = (0 - B);
C = (0 - C);
}
normalize_i64_i64_i64_ptr_i64_ptr_i64_ptr_i64(A, B, C, box_a, box_b, box_c);
aa = (((unsigned)(0) < 1) ? box_a[0] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(0), 1), flow_fault_handler("array index out of bounds"), box_a[0]));
bb = (((unsigned)(0) < 1) ? box_b[0] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(0), 1), flow_fault_handler("array index out of bounds"), box_b[0]));
cc = (((unsigned)(0) < 1) ? box_c[0] : (fprintf(stderr, "array index %d out of bounds (size %d)\n", (int)(0), 1), flow_fault_handler("array index out of bounds"), box_c[0]));
nn = m;
sign = (0 - sign);
}
return res;
}
int64_t P_func_i64_i64_i64_i64_i64_i64_i64(int64_t k, int64_t L, int64_t U, int64_t len_side, int64_t total_pairs, int64_t U2m1, int64_t Lm1) {
if (k == 1) {
return FLOW_CHECKED_DIV(((len_side * (len_side - 1))), (2));
}
int64_t r = isqrt_i64(k);
if ((r * r) == k) {
int64_t s = r;
int64_t b = FLOW_CHECKED_DIV((U), (s));
if (b < L) {
return total_pairs;
}
if (b > U) {
b = U;
}
int64_t cnt = ((b - L) + 1);
int64_t sum_x = FLOW_CHECKED_DIV((((L + b) * cnt)), (2));
int64_t partial = ((s * sum_x) - (L * cnt));
int64_t full_cnt = (U - b);
return (partial + (full_cnt * len_side));
}
int64_t b2 = isqrt_i64(FLOW_CHECKED_DIV((U2m1), (k)));
if (b2 < L) {
return total_pairs;
}
if (b2 > U) {
b2 = U;
}
int64_t cnt2 = ((b2 - L) + 1);
int64_t sum_floor = (sum_floor_mul_surd_i64_i64_i64_i64_i64(0, 1, 1, k, b2) - sum_floor_mul_surd_i64_i64_i64_i64_i64(0, 1, 1, k, Lm1));
int64_t partial2 = (sum_floor - ((L - 1) * cnt2));
int64_t full_cnt2 = (U - b2);
return (partial2 + (full_cnt2 * len_side));
}
int64_t solve_i64_i64(int64_t M, int64_t N) {
int64_t L = (M + 1);
int64_t U = N;
if (L > U) {
return 0;
}
int64_t len_side = ((U - L) + 1);
int64_t total_pairs = (len_side * len_side);
int64_t k_max = FLOW_CHECKED_DIV(((U * U)), ((L * L)));
int64_t last_odd = k_max;
if ((k_max & 1) == 0) {
last_odd = (k_max - 1);
}
int64_t U2m1 = ((U * U) - 1);
int64_t Lm1 = (L - 1);
int64_t ans = 0;
int64_t n = 1;
while (n <= last_odd) {
ans = ((ans + P_func_i64_i64_i64_i64_i64_i64_i64((n + 1), L, U, len_side, total_pairs, U2m1, Lm1)) - P_func_i64_i64_i64_i64_i64_i64_i64(n, L, U, len_side, total_pairs, U2m1, Lm1));
n = (n + 2);
}
return ans;
}
int32_t main(void) {
int64_t ans = solve_i64_i64(2000000, 1000000000);
printf("%lld\n", ans);
return 0;
}