# Project Euler 332
# Sum of minimal spherical triangle areas A(r) for r=1..50.
import euler.nt { isqrt }
extern {
function calloc(n: i64, size: i64) -> ptr<void>
function free(p: ptr<void>) -> void
function atan2(y: f64, x: f64) -> f64
function tan(x: f64) -> f64
function floor(x: f64) -> f64
}
function gcd_abs(a0: i64, b0: i64) -> i64 {
let mut a: i64 = a0
let mut b: i64 = b0
if a < 0 { a = -a }
if b < 0 { b = -b }
while b != 0 {
let t: i64 = a % b
a = b
b = t
}
return a
}
function min_area(r: i64) -> f64 {
let rr: i64 = r * r
let xs: ptr<i16> = calloc(30000, 2)
let ys: ptr<i16> = calloc(30000, 2)
let zs: ptr<i16> = calloc(30000, 2)
let mut n: i64 = 0
let mut x: i64 = -r
while x <= r {
let x2: i64 = x * x
let mut y: i64 = -r
while y <= r {
let z2: i64 = rr - x2 - y * y
if z2 >= 0 {
let z: i64 = isqrt(z2)
if z * z == z2 {
xs[n] = x as i16; ys[n] = y as i16; zs[n] = z as i16; n = n + 1
if z != 0 {
xs[n] = x as i16; ys[n] = y as i16; zs[n] = (-z) as i16; n = n + 1
}
}
}
y = y + 1
}
x = x + 1
}
if n < 3 { free(xs); free(ys); free(zs); return 0.0 }
let dots: ptr<i16> = calloc(n * n, 2)
let mut i: i64 = 0
while i < n {
let xi: i64 = xs[i] as i64
let yi: i64 = ys[i] as i64
let zi: i64 = zs[i] as i64
let mut j: i64 = i
while j < n {
let v: i64 = xi * (xs[j] as i64) + yi * (ys[j] as i64) + zi * (zs[j] as i64)
dots[i * n + j] = v as i16
dots[j * n + i] = v as i16
j = j + 1
}
i = i + 1
}
let r2: f64 = (rr) as f64
let r3: f64 = r2 * (r as f64)
let inv_r2: f64 = 1.0 / r2
let inv_r3: f64 = 1.0 / r3
let mut best: f64 = 1.5707963267948966 * r2
let mut best_unit: f64 = best / r2
let mut D_max: i64 = floor(4.0 * r3 * tan(best_unit / 2.0) + 0.000000000001) as i64
if D_max < 1 { D_max = 1 }
i = 0
while i < n - 2 {
let xi: i64 = xs[i] as i64
let yi: i64 = ys[i] as i64
let zi: i64 = zs[i] as i64
let mut j: i64 = i + 1
while j < n - 1 {
let xj: i64 = xs[j] as i64
let yj: i64 = ys[j] as i64
let zj: i64 = zs[j] as i64
let cx: i64 = yi * zj - zi * yj
let cy: i64 = zi * xj - xi * zj
let cz: i64 = xi * yj - yi * xj
if cx != 0 || cy != 0 || cz != 0 {
let g: i64 = gcd_abs(cx, gcd_abs(cy, cz))
if g != 0 && g <= D_max {
let dot_ij: i64 = dots[i * n + j] as i64
let mut k: i64 = j + 1
while k < n {
let mut det: i64 = cx * (xs[k] as i64) + cy * (ys[k] as i64) + cz * (zs[k] as i64)
if det != 0 {
if det < 0 { det = -det }
if det <= D_max {
let denom: f64 = 1.0 + ((dot_ij + (dots[j * n + k] as i64) + (dots[i * n + k] as i64)) as f64) * inv_r2
let area: f64 = 2.0 * atan2((det as f64) * inv_r3, denom) * r2
if area < best {
best = area
best_unit = best / r2
D_max = floor(4.0 * r3 * tan(best_unit / 2.0) + 0.000000000001) as i64
if D_max < 1 { D_max = 1 }
}
}
}
k = k + 1
}
}
}
j = j + 1
}
i = i + 1
}
free(xs); free(ys); free(zs); free(dots)
return best
}
function main() -> i32 {
let mut sum: f64 = 0.0
let mut r: i64 = 1
while r <= 50 {
sum = sum + min_area(r)
r = r + 1
}
printf("%.6f\n", sum)
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);
double atan2(double y, double x);
int64_t gcd_abs_i64_i64(int64_t a0, int64_t b0);
double min_area_i64(int64_t r);
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 gcd_abs_i64_i64(int64_t a0, int64_t b0) {
int64_t a = a0;
int64_t b = b0;
if (a < 0) {
a = (-a);
}
if (b < 0) {
b = (-b);
}
while (b != 0) {
int64_t t = FLOW_CHECKED_MOD((a), (b));
a = b;
b = t;
}
return a;
}
double min_area_i64(int64_t r) {
int64_t rr = (r * r);
int16_t* xs = (int16_t*)(calloc(30000, 2));
int16_t* ys = (int16_t*)(calloc(30000, 2));
int16_t* zs = (int16_t*)(calloc(30000, 2));
int64_t n = 0;
int64_t x = (-r);
while (x <= r) {
int64_t x2 = (x * x);
int64_t y = (-r);
while (y <= r) {
int64_t z2 = ((rr - x2) - (y * y));
if (z2 >= 0) {
int64_t z = isqrt_i64(z2);
if ((z * z) == z2) {
xs[n] = ((int16_t)(x));
ys[n] = ((int16_t)(y));
zs[n] = ((int16_t)(z));
n = (n + 1);
if (z != 0) {
xs[n] = ((int16_t)(x));
ys[n] = ((int16_t)(y));
zs[n] = ((int16_t)((-z)));
n = (n + 1);
}
}
}
y = (y + 1);
}
x = (x + 1);
}
if (n < 3) {
free(xs);
free(ys);
free(zs);
return 0.0;
}
int16_t* dots = (int16_t*)(calloc((n * n), 2));
int64_t i = 0;
while (i < n) {
int64_t xi = ((int64_t)(xs[i]));
int64_t yi = ((int64_t)(ys[i]));
int64_t zi = ((int64_t)(zs[i]));
int64_t j = i;
while (j < n) {
int64_t v = (((xi * ((int64_t)(xs[j]))) + (yi * ((int64_t)(ys[j])))) + (zi * ((int64_t)(zs[j]))));
dots[((i * n) + j)] = ((int16_t)(v));
dots[((j * n) + i)] = ((int16_t)(v));
j = (j + 1);
}
i = (i + 1);
}
double r2 = ((double)(rr));
double r3 = (r2 * ((double)(r)));
double inv_r2 = (1.0 / r2);
double inv_r3 = (1.0 / r3);
double best = (1.5707963267948966 * r2);
double best_unit = (best / r2);
int64_t D_max = ((int64_t)(floor((((4.0 * r3) * tan((best_unit / 2.0))) + 0.000000000001))));
if (D_max < 1) {
D_max = 1;
}
i = 0;
while (i < (n - 2)) {
int64_t xi = ((int64_t)(xs[i]));
int64_t yi = ((int64_t)(ys[i]));
int64_t zi = ((int64_t)(zs[i]));
int64_t j = (i + 1);
while (j < (n - 1)) {
int64_t xj = ((int64_t)(xs[j]));
int64_t yj = ((int64_t)(ys[j]));
int64_t zj = ((int64_t)(zs[j]));
int64_t cx = ((yi * zj) - (zi * yj));
int64_t cy = ((zi * xj) - (xi * zj));
int64_t cz = ((xi * yj) - (yi * xj));
if (((cx != 0 || cy != 0) || cz != 0)) {
int64_t g = gcd_abs_i64_i64(cx, gcd_abs_i64_i64(cy, cz));
if ((g != 0 && g <= D_max)) {
int64_t dot_ij = ((int64_t)(dots[((i * n) + j)]));
int64_t k = (j + 1);
while (k < n) {
int64_t det = (((cx * ((int64_t)(xs[k]))) + (cy * ((int64_t)(ys[k])))) + (cz * ((int64_t)(zs[k]))));
if (det != 0) {
if (det < 0) {
det = (-det);
}
if (det <= D_max) {
double denom = (1.0 + (((double)(((dot_ij + ((int64_t)(dots[((j * n) + k)]))) + ((int64_t)(dots[((i * n) + k)]))))) * inv_r2));
double area = ((2.0 * atan2((((double)(det)) * inv_r3), denom)) * r2);
if (area < best) {
best = area;
best_unit = (best / r2);
D_max = ((int64_t)(floor((((4.0 * r3) * tan((best_unit / 2.0))) + 0.000000000001))));
if (D_max < 1) {
D_max = 1;
}
}
}
}
k = (k + 1);
}
}
}
j = (j + 1);
}
i = (i + 1);
}
free(xs);
free(ys);
free(zs);
free(dots);
return best;
}
int32_t main(void) {
double sum = 0.0;
int64_t r = 1;
while (r <= 50) {
sum = (sum + min_area_i64(r));
r = (r + 1);
}
printf("%.6f\n", sum);
return 0;
}