Asymmetric Diophantine Equation: 16x^2 + y^4 = z^2 with x, y, z positive integers and gcd(x, y, z) = 1. S(N) sums x + y + z over all solutions with 1 <= x, y, z <= N. Find S(10^16) mod 10^9. Rewrite as (4x)^2 + (y^2)^2 = z^2, a Pythagorean triple. y odd: the triple is primitive, so y^2 = m^2 - n^2, 4x = 2mn, z = m^2 + n^2. Then m - n = s^2, m + n = t^2 for odd coprime s < t, giving x = (t^4 - s^4)/8, y = st, z = (t^4 + s^4)/2. y even: x is odd, z = 4w, and (x, (y/2)^2, w) is primitive: x = m^2 - n^2, (y/2)^2 = 2mn, w = m^2 + n^2. Since 2mn is a square with gcd(m, n) = 1, {m, n} = {2s^2, t^2} with t odd and gcd(s, t) = 1, so y = 4st and, by which of 2s^2, t^2 is larger, x = 4s^4 - t^4, z = 16s^4 + 4t^4 (2s^2 > t^2), or x = t^4 - 4s^4, z = 4t^4 + 16s^4 (t^2 > 2s^2). Each primitive triple arises once, so the three families enumerate every solution exactly once. z <= 10^16 bounds s, t near 12000; everything fits in i64 and the running sum is reduced mod 10^9.
# Project Euler 764
# Asymmetric Diophantine Equation: 16x^2 + y^4 = z^2 with x, y, z positive
# integers and gcd(x, y, z) = 1. S(N) sums x + y + z over all solutions
# with 1 <= x, y, z <= N. Find S(10^16) mod 10^9.
#
# Rewrite as (4x)^2 + (y^2)^2 = z^2, a Pythagorean triple.
# y odd: the triple is primitive, so y^2 = m^2 - n^2, 4x = 2mn,
# z = m^2 + n^2. Then m - n = s^2, m + n = t^2 for odd coprime
# s < t, giving x = (t^4 - s^4)/8, y = st, z = (t^4 + s^4)/2.
# y even: x is odd, z = 4w, and (x, (y/2)^2, w) is primitive:
# x = m^2 - n^2, (y/2)^2 = 2mn, w = m^2 + n^2. Since 2mn is a
# square with gcd(m, n) = 1, {m, n} = {2s^2, t^2} with t odd and
# gcd(s, t) = 1, so y = 4st and, by which of 2s^2, t^2 is larger,
# x = 4s^4 - t^4, z = 16s^4 + 4t^4 (2s^2 > t^2), or
# x = t^4 - 4s^4, z = 4t^4 + 16s^4 (t^2 > 2s^2).
# Each primitive triple arises once, so the three families enumerate every
# solution exactly once. z <= 10^16 bounds s, t near 12000; everything
# fits in i64 and the running sum is reduced mod 10^9.
const N: i64 = 10000000000000000
const MOD: i64 = 1000000000
function gcd(a: i64, b: i64) -> i64 {
let mut x: i64 = a
let mut y: i64 = b
while y != 0 {
let r: i64 = x % y
x = y
y = r
}
return x
}
function main() -> i32 {
let mut total: i64 = 0
# Family 1: y odd. Odd coprime s < t.
# x = (t^4 - s^4)/8, y = s t, z = (t^4 + s^4)/2.
let mut t: i64 = 3
while (t * t * t * t + 1) / 2 <= N {
let t4: i64 = t * t * t * t
let mut s: i64 = 1
while s < t {
let s4: i64 = s * s * s * s
let z: i64 = (t4 + s4) / 2
if z > N {
s = t
} else {
if gcd(s, t) == 1 {
let x: i64 = (t4 - s4) / 8
let y: i64 = s * t
if x >= 1 && x <= N && y <= N {
total = (total + (x + y + z) % MOD) % MOD
}
}
s = s + 2
}
}
t = t + 2
}
# Family 2: y even, m = 2s^2 > n = t^2, t odd, gcd(s, t) = 1.
# x = 4s^4 - t^4, y = 4 s t, z = 16s^4 + 4t^4.
let mut s2: i64 = 1
while 16 * s2 * s2 * s2 * s2 + 4 <= N {
let s4: i64 = s2 * s2 * s2 * s2
let mut tt: i64 = 1
while tt * tt < 2 * s2 * s2 {
let t4: i64 = tt * tt * tt * tt
let z: i64 = 16 * s4 + 4 * t4
if z > N {
tt = 2 * s2 * s2
} else {
if gcd(s2, tt) == 1 {
let x: i64 = 4 * s4 - t4
let y: i64 = 4 * s2 * tt
if x >= 1 && x <= N && y <= N {
total = (total + (x + y + z) % MOD) % MOD
}
}
tt = tt + 2
}
}
s2 = s2 + 1
}
# Family 3: y even, m = t^2 > n = 2s^2, t odd, gcd(s, t) = 1.
# x = t^4 - 4s^4, y = 4 s t, z = 4t^4 + 16s^4.
let mut t3: i64 = 1
while 4 * t3 * t3 * t3 * t3 + 16 <= N {
let t4: i64 = t3 * t3 * t3 * t3
let mut ss: i64 = 1
while 2 * ss * ss < t3 * t3 {
let s4b: i64 = ss * ss * ss * ss
let z: i64 = 4 * t4 + 16 * s4b
if z > N {
ss = t3 * t3
} else {
if gcd(ss, t3) == 1 {
let x: i64 = t4 - 4 * s4b
let y: i64 = 4 * ss * t3
if x >= 1 && x <= N && y <= N {
total = (total + (x + y + z) % MOD) % MOD
}
}
ss = ss + 1
}
}
t3 = t3 + 2
}
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; }
int64_t gcd_i64_i64(int64_t a, int64_t b);
int32_t main(void);
static const int64_t N = 10000000000000000;
static const int64_t MOD = 1000000000;
int64_t gcd_i64_i64(int64_t a, int64_t b) {
int64_t x = a;
int64_t y = b;
while (y != 0) {
int64_t r = FLOW_CHECKED_MOD((x), (y));
x = y;
y = r;
}
return x;
}
int32_t main(void) {
int64_t total = 0;
int64_t t = 3;
while (FLOW_CHECKED_DIV((((((t * t) * t) * t) + 1)), (2)) <= N) {
int64_t t4 = (((t * t) * t) * t);
int64_t s = 1;
while (s < t) {
int64_t s4 = (((s * s) * s) * s);
int64_t z = FLOW_CHECKED_DIV(((t4 + s4)), (2));
if (z > N) {
s = t;
} else {
if (gcd_i64_i64(s, t) == 1) {
int64_t x = FLOW_CHECKED_DIV(((t4 - s4)), (8));
int64_t y = (s * t);
if (((x >= 1 && x <= N) && y <= N)) {
total = FLOW_CHECKED_MOD(((total + FLOW_CHECKED_MOD((((x + y) + z)), (MOD)))), (MOD));
}
}
s = (s + 2);
}
}
t = (t + 2);
}
int64_t s2 = 1;
while ((((((16 * s2) * s2) * s2) * s2) + 4) <= N) {
int64_t s4 = (((s2 * s2) * s2) * s2);
int64_t tt = 1;
while ((tt * tt) < ((2 * s2) * s2)) {
int64_t t4 = (((tt * tt) * tt) * tt);
int64_t z = ((16 * s4) + (4 * t4));
if (z > N) {
tt = ((2 * s2) * s2);
} else {
if (gcd_i64_i64(s2, tt) == 1) {
int64_t x = ((4 * s4) - t4);
int64_t y = ((4 * s2) * tt);
if (((x >= 1 && x <= N) && y <= N)) {
total = FLOW_CHECKED_MOD(((total + FLOW_CHECKED_MOD((((x + y) + z)), (MOD)))), (MOD));
}
}
tt = (tt + 2);
}
}
s2 = (s2 + 1);
}
int64_t t3 = 1;
while ((((((4 * t3) * t3) * t3) * t3) + 16) <= N) {
int64_t t4 = (((t3 * t3) * t3) * t3);
int64_t ss = 1;
while (((2 * ss) * ss) < (t3 * t3)) {
int64_t s4b = (((ss * ss) * ss) * ss);
int64_t z = ((4 * t4) + (16 * s4b));
if (z > N) {
ss = (t3 * t3);
} else {
if (gcd_i64_i64(ss, t3) == 1) {
int64_t x = (t4 - (4 * s4b));
int64_t y = ((4 * ss) * t3);
if (((x >= 1 && x <= N) && y <= N)) {
total = FLOW_CHECKED_MOD(((total + FLOW_CHECKED_MOD((((x + y) + z)), (MOD)))), (MOD));
}
}
ss = (ss + 1);
}
}
t3 = (t3 + 2);
}
printf("%lld\n", total);
return 0;
}