#include <chrono>
#include <concepts>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <stdexcept>
#include "SafeInt.hpp"

/* Martin Fuller 2025-03-16
Proof that a(n) > 2*10^(n-1) + ((2/3)*10^(n-1))^0.5:

Let m = 10^(n-1). We know that a,b,c,d >= m and that e >= 2m.
Rewrite the quintuple as: (m+p)^2 + (m+q)^2 + (m+r)^2 + (m+s)^2 = (2m+t)^2
Expand then cancel the term in m^2: 2m(p + q + r + s) + (p^2 + q^2 + r^2 + s^2) = 2m(2t) + t^2

Lemma: 2t > p+q+r+s
Suppose to the contrary that p+q+r+s >= 2t:
Then the mean of p,q,r,s is at least t/2, and root-mean-square is the same or larger.
Equality is not possible since p,q,r,s cannot all be equal to t/2.
Therefore p^2 + q^2 + r^2 + s^2 > t^2, which contradicts the original quintuple equation.

Let f = 2t - (p+q+r+s) >= 1.  Then p^2 + q^2 + r^2 + s^2 = 2mf + t^2.
Substitute t^2 > ((p+q+r+s)/2)^2 > (p^2 + q^2 + r^2 + s^2) / 4
=> p^2 + q^2 + r^2 + s^2 > (8/3)mf

Create a lower bound for p+q+r+s by setting p,q,r=0 and s^2=(8/3)mf
t > (p+q+r+s)/2 > ((2/3)mf)^0.5 >= ((2/3)m)^0.5
a(n) = e = 2m+t > 2*10^(n-1) + ((2/3)*10^(n-1))^0.5
//

Background to the algorithm:

Given m,t,f,p,q we can calculate r,s as follows:
  Let y = r+s = 2t-f-p-q and z = r^2+s^2 = 2mf+t^2-p^2-q^2
  z = r^2+(y-r)^2
  => r,s = y/2 -+ sqrt(2z - y^2)/2

Given m,t,f,p we can calculate bounds on q as follows:
  Let w = q+r+s = 2t-f-p and x = q^2+r^2+s^2 = 2mf+t^2-p^2
  Note: y = w-q and z = x-q^2.

  Require q < w/3 because q+r+s=w and q<r<s.

  Require 2z-y^2 >= 0 for r,s to be real.
  -3q^2 + 2wq + 2x - w^2 >= 0
  => q >= w/3 - sqrt(6x - 2w^2)/3

  Also require q < r = y/2-sqrt(2z-y^2)/2
  y-2q > sqrt(2z-y^2)
  (w-3q)^2 > 2(x-q^2) - (w-q)^2
  12q^2 - 8wq - 2x + 2w^2 > 0
  => q < w/3 - sqrt(6x - 2w^2)/6

Given m,t,f we can calculate bounds on p as follows:
  Let u = p+q+r+s = 2t-f and v = p^2+q^2+r^2+s^2 = 2mf+t^2

  Require 0 <= p < u/4 because p+q+r+s=u and p<q<r<s

  Also require p < q < w/3 - sqrt(6x - 2w^2)/6
  (6x - 2w^2)/36 < (w/3 - p)^2
  12p^2 - 6up + u^2 - v > 0
  => p < u/4 - sqrt(12v - 3u^2)/12

Given m,t we can calculate bounds on f as follows:
  t > ((2/3)mf)^0.5
  => f < (3/2)t^2 / m

  Also require 0 <= u/4 - sqrt(12v - 3u^2)/12
  Notice that u is a decreasing function of f, and v is an increasing function of f.
  => If the inequality is not satisfied for some f, then it is also not satisfied for larger f.

Lemma: e cannot be a multiple of 4.
Proof: The squares mod 8 are {0,1,4}. At least one of a,b,c,d is odd, because they have no common factor. Therefore a^2+b^2+c^2+d^2 cannot be 0 mod 8.
Similar: If e is even then a,b,c,d must all be odd.

The algorithm has 4 nested loops: t in ascending order from the lower bound, then f, then p, then q. The inner loop calculates r,s and checks if they are positive integers and that a,b,c,d are in ascending order with no common factor.
The program is implemented using 64-bit unsigned integer arithmetic. It uses SafeInt to guard against overflow/underflow.
It is relatively easy to use 64-bit arithmetic up to n=19 based on p,q,r,s,t because they are are significantly smaller than a,b,c,d,e.
*/

typedef uint_fast64_t a380729int;
typedef SafeInt<a380729int> a380729safe;

template <typename T>
std::ostream& operator<<(std::ostream &os, const SafeInt<T>& t) {
  return os << static_cast<T>(t);
}

template <typename T>
T constexpr ipow(T x, std::unsigned_integral auto n) {
  if (n == 0) return 1;
  T y = 1;
  while (n > 1) {
    if (n & 1) y *= x;
    x *= x;
    n >>= 1;
  }
  return x * y;
}

void test_ipow() {
  std::cerr << "ipow\n";
  for (unsigned n = 0; n < 64; n++) {
    if (ipow(0u, n) != !n) { std::cerr << "ipow(0u," << n << ")==" << ipow(0u,n) << '\n'; return; }
  }
  for (unsigned n = 0; n < 64; n++) {
    if (ipow(1u, n) != 1u) { std::cerr << "ipow(1u," << n << ")==" << ipow(1u,n) << '\n'; return; }
  }
  for (uint64_t x = 2; x <= 64; x++) {
    uint64_t y = 1, limit = std::numeric_limits<uint64_t>::max()/x;
    for (unsigned n = 0; y <= limit; n++, y*=x) {
      if (ipow(x,n) != y) { std::cerr << "ipow(" << x << ',' << n << ")==" << ipow(x,n) << '\n'; return; }
    }
  }
}

template <typename U>
U constexpr isqrt(U x) {
  U r = static_cast<U>(std::sqrt(static_cast<double>(x))), s = r*r;
  if (s <= x && s >= r) {
    while (s + 2*r < x) {
      s += 2*r + 1;
      r++;
    }
  } else {
    do {
      r--;
      s -= 2*r+1;
    } while (s > x);
  }
  return r;
}

void test_isqrt() {
  std::cerr << "isqrt\n";
  auto start = std::chrono::steady_clock::now();
  for (uint64_t x = 0, y = 0; x < uint64_t{1} << 32; y+=2*x+1, x++) {
    if (isqrt(y) != x) { std::cerr << "isqrt(" << y << ")==" << isqrt(y) << '\n'; return; }
    if (isqrt(y+2*x) != x) { std::cerr << "isqrt(" << (y+2*x) << ")==" << isqrt(y+2*x) << '\n'; return; }
  }
  auto duration = std::chrono::steady_clock::now() - start;
  std::cerr << "isqrt " << std::chrono::duration_cast<std::chrono::seconds>(duration) << '\n';
}

auto constexpr divUp(a380729safe n, a380729safe d) {
  return (n / d) + ((n % d) ? 1 : 0);
}

template <typename U>
U constexpr isqrtUp(U x) {
  if (!x) return 0;
  return isqrt(x-U{1})+U{1};
}

a380729int gcd(a380729int a, a380729int b, a380729int c, a380729int d) {
  return std::gcd(std::gcd(std::gcd(a, b), c), d);
}

a380729int a380729(unsigned n) {
  if (!n) throw std::domain_error("n == 0");
  const a380729safe m = ipow(a380729safe{10}, n-1);
  for (a380729safe t = isqrtUp(m - m/3); true; t++) {
    if (((2*m+t) & 3) == 0) continue;
    std::cerr << static_cast<a380729int>(t) << ',';
    a380729safe fMax = t*t*3/(2*m); 
    for (unsigned f = 1; f <= fMax; f++) {
      const a380729safe u = 2*t-f, v = 2*m*f + t*t;
      const a380729safe pMin = 0;
      const a380729safe pTemp = isqrt(3*v - 3*u*u/4) * 2; /* rounded down from isqrt(12*v - 3*u*u) */
      if (pTemp > 3*u) break;
      const a380729safe pMax = (3*u - pTemp) / 12;
      for (a380729safe p = pMin; p <= pMax; p++) {
        if (!(t&1) && !((m+p)&1)) continue;
        const a380729safe w = u-p, x = v-p*p;
        if (6*x < 2*w*w) continue;
        const a380729safe qMin = std::max(p+1, divUp(w - std::min(w, isqrt(6*x - 2*w*w)), 3u));
        const a380729safe qMax = (w*2 - std::min(w*2, isqrtUp(6*x - 2*w*w))) / 6;
        for (a380729safe q = qMin; q <= qMax; q++) {
          const a380729safe y = w-q, z = x-q*q;
          const a380729safe discr = 2*z-y*y;
          const a380729safe sqrt = isqrt(discr);
          if (sqrt*sqrt == discr) {
            a380729safe r = y - sqrt;
            if (r & 1) continue;
            r /= 2;
            a380729safe s = y - r;
            if (q < r && r < s && gcd(m+p,m+q,m+r,m+s) == 1) {
              std::cerr << '\n';
              std::cout << '[' << (m+p) << ", " << (m+q) << ", " << (m+r) << ", " << (m+s) << ", " << (2*m+t) << "]\n";
              return 2*m+t;
            }
          }
        }
      }
    }
  }
}

int main() {
//  test_ipow();
//  test_isqrt();
  for (unsigned n = 1; n <= 19; n++) {
    a380729(n);
  }
}