pythagorean-triple.c

#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

typedef uint64_t n;

n sqrt_n(n x) {
  if (x <= 1) return x;

  n x0 = x / 2;
  n x1 = (x0 + x / x0) / 2;

  while (x1 < x0) {
    x0 = x1;
    x1 = (x0 + x / x0) / 2;
  }
  return x0;
}

bool check(n x, n y, n *z) {
  n x2 = x * x;
  if (x2 < x) return false;
  n y2 = y * y;
  if (y2 < y) return false;
  n sum = x2 + y2;
  if (sum < x || sum < y) return false;

  *z = sqrt_n(sum);

  return *z * *z == sum;
}

bool is_prime(n x) {
  if (x <= 1) return false;
  if (x % 2 == 0 && x > 2) return false;
  for (n i = 3; i <= sqrt_n(x); i += 2) {
    if (x % i == 0) return false;
  }
  return true;
}

int main(int argc, char **argv) {
  if (argc != 2) return 1;
  n max;
  if (sscanf(argv[1], "%llu", &max) != 1) return 1;
  n z;

  for (n x = 1; x < max; ++x) {
    for (n y = x + 1; y < max; ++y) {
      if (check(x, y, &z) && (is_prime(x) || is_prime(y) || is_prime(z))) {
        printf("%llu %llu %llu\n", x, y, z);
      }
    }
  }
  return 0;
}