# Written in 2025 by Kirill A. Korinsky <kirill@korins.ky>
#
# To the extent possible under law, the author(s) have dedicated all
# copyright and related and neighboring rights to this software to the
# public domain worldwide.
#
# This software is distributed without any warranty.

import itertools
import math
from typing import Iterator, Optional, Tuple


class KunerthSolver:
    """
    Deterministic solver for quadratic Diophantine equations
    y^2 = ax^2 + bx + c.
    Adolf Kunerth's 1877-1880 algorithms.

    1. Modular Case (a == 0):
       - Shift strategy c' = c + kb iterates until Jacobi symbol
         confirms ±b is quadratic residue modulo c'.
       - Continued fraction convergents of sqrt(c') generate parameters
         t minimizing |t^2 - c'|.

    2. General Case (a != 0):
       - Outer loop: Wavefront enumeration of (alpha, beta) for square
         injection.
       - Inner loop: Continued fraction convergents of sqrt(a) solve
         v^2 - aw^2 = N where N divides Stem Number.
    """

    def __init__(self, a: int, b: int, c: int):
        self.a = a
        self.b = b
        self.c = c

    @staticmethod
    def _jacobi(a: int, n: int) -> int:
        if n <= 0 or n % 2 == 0:
            return 0

        a %= n
        result = 1
        while a != 0:
            while a % 2 == 0:
                a //= 2
                n_mod_8 = n % 8
                if n_mod_8 in (3, 5):
                    result = -result
            a, n = n, a
            if a % 4 == 3 and n % 4 == 3:
                result = -result
            a %= n
        if n == 1:
            return result
        return 0

    def _check_hauptformel(self, seed: Tuple[int, int, int],
                           param: Tuple[int, int]) -> Optional[Tuple[int, int]]:
        m, n, r = seed
        v, w = param

        # Denominator: v^2 - a*w^2 (Formula X denominator)
        denom = v * v - self.a * w * w
        if denom == 0:
            return None

        # Constant K = 2am + bn (Derived from Formula VIII numerator terms)
        K = 2 * self.a * m + self.b * n

        # Stem Number Check (Kunerth 1880)
        # S = (2am + bn)^2 - 4ar^2 = K^2 - 4ar^2
        # The denominator must divide the Stem Number.
        stem = K * K - 4 * self.a * r * r
        if stem % denom != 0:
            return None

        # Numerator term: w * (2rv + Kw)
        # Checks if the rational expression yields an integer
        num_term = w * (2 * r * v + K * w)

        if num_term % denom != 0:
            return None
        term = num_term // denom

        # Calculate x: nx = m + term
        rhs = m + term
        if rhs % n != 0:
            return None
        x = rhs // n

        # Recover y via Formula IX: y = p(x - m/n) - r/n
        # Simplified to integer arithmetic: y = (v(nx - m) - rw) / nw
        y_num = v * term - r * w
        y_den = n * w

        if y_den == 0:
            return None
        if y_num % y_den != 0:
            return None
        y = y_num // y_den

        # Verify
        if y * y == self.a * x * x + self.b * x + self.c:
            return x, y
        if (-y)**2 == self.a * x * x + self.b * x + self.c:
            return x, -y

        return None

    def _get_modular_seed_components(self, t: int) -> Tuple[int, int, int]:
        num = t * t - self.c
        den = self.b
        g = math.gcd(num, den)
        m, n = num // g, den // g
        if n < 0:
            m, n = -m, -n
        r = n * t
        return (m, n, r)

    @staticmethod
    def _continued_fraction_gen(r: int) -> Iterator[Tuple[int, int]]:
        if r <= 0:
            return
        a0 = int(math.isqrt(r))
        if a0 * a0 == r:
            yield (a0, 1)
            return

        h1, h2 = 1, 0
        k1, k2 = 0, 1

        m = 0
        d = 1
        a_n = a0

        stop_coeff = 2 * a0

        while True:
            h = a_n * h1 + h2
            k = a_n * k1 + k2

            yield (h, k)

            h2, h1 = h1, h
            k2, k1 = k1, k

            m = d * a_n - m
            d = (r - m * m) // d

            if d == 0:
                break

            a_n = (a0 + m) // d

            # The period of sqrt(r) always ends with 2*a0.
            if a_n == stop_coeff:
                h_final = a_n * h1 + h2
                k_final = a_n * k1 + k2
                yield (h_final, k_final)
                return

    def _cf_shifted_strategy(self) -> Iterator[Tuple[int, int]]:
        # We need c + k*b > 0  =>  k*b > -c  =>  k > -c / b
        if self.c < 0:
            start_k = (abs(self.c) // self.b) + 1
        else:
            start_k = 0

        v_variants = [1, -1, 2, -2]

        # We have no idea about b strucutre, so, jacobi may lie
        k_lim = abs(self.b)
        t_lim = abs(self.b)

        for k in itertools.count(start_k):
            if k > k_lim:
                break

            c_prime = self.c + k * self.b

            if c_prime <= 0:
                continue

            allow_pos = False
            allow_neg = False

            jacobi_mod = c_prime
            while jacobi_mod % 2 == 0:
                jacobi_mod //= 2

            if jacobi_mod == 1:
                allow_pos = True
                allow_neg = True
            else:
                if self._jacobi(self.b, jacobi_mod) != -1:
                    allow_pos = True
                if self._jacobi(-self.b, jacobi_mod) != -1:
                    allow_neg = True

            if not allow_pos and not allow_neg:
                continue

            cf_gen = self._continued_fraction_gen(c_prime)

            for t, _ in cf_gen:
                if abs(t) > t_lim:
                    break

                diff = t * t - c_prime

                if diff > 0 and not allow_pos:
                    continue
                if diff < 0 and not allow_neg:
                    continue

                seed = self._get_modular_seed_components(t)

                for v in v_variants:
                    res = self._check_hauptformel(seed, (v, 1))
                    if res:
                        yield res

    def _wavefront_pairs(self) -> Iterator[Tuple[int, int]]:
        root_a = int(math.isqrt(max(0, self.a)))
        dist = 0
        while True:
            # Generate pairs (alpha, beta) in a wavefront
            for i in range(-dist, dist + 1):
                yield (root_a + i, dist)
                if dist > 0:
                    yield (root_a + i, -dist)

            for j in range(-dist + 1, dist):
                yield (root_a + dist, j)
                yield (root_a - dist, j)

            dist += 1

    def _wavefront_seeds(self) -> Iterator[Tuple[int, int, int]]:
        for alpha, beta in self._wavefront_pairs():
            # Coefficients of the remainder polynomial R(x) (Kunerth 1878, Eq [2])
            A_p = self.a - alpha**2
            B_p = self.b - 2 * alpha * beta
            C_p = self.c - beta**2

            # Case 1: Linear Remainder (Equation matches format y^2 = (ax+b)^2 + linear)
            if A_p == 0:
                if B_p == 0:
                    continue
                g = math.gcd(C_p, B_p)
                m, n = -C_p // g, B_p // g
                if n < 0:
                    m, n = -m, -n
                r = alpha * m + beta * n
                if n != 0 and r * r == self.a * m * m + self.b * m * n + self.c * n * n:
                    yield m, n, r
                continue

            # Case 2: Quadratic Remainder (Discriminant condition, Kunerth 1878, Eq [3])
            delta = B_p**2 - 4 * A_p * C_p
            if delta < 0:
                continue
            s_delta = math.isqrt(delta)
            if s_delta**2 != delta:
                continue

            for sign in (1, -1):
                num = -B_p + sign * s_delta
                den = 2 * A_p
                if den == 0:
                    continue
                g = math.gcd(num, den)
                m, n = num // g, den // g
                if n < 0:
                    m, n = -m, -n
                r = alpha * m + beta * n
                if r * r == self.a * m * m + self.b * m * n + self.c * n * n:
                    yield m, n, r

    def solve(self) -> Optional[Tuple[int, int]]:
        # Early Exit Criteria (Kunerth 1880, Achterreste)
        if self.a == 0:
            if self.b % 4 == 0 and self.c % 4 in (2, 3):
                return None

            if self.b % 8 == 0 and self.c % 8 == 5:
                return None

        if self.a == 0:
            gen_cf = self._cf_shifted_strategy()

            for res in gen_cf:
                return res

        else:
            seed_gen = self._wavefront_seeds()

            for seed in seed_gen:
                cf_gen = self._continued_fraction_gen(self.a)

                for param in cf_gen:
                    res = self._check_hauptformel(seed, param)
                    if res:
                        return res

        return None


if __name__ == "__main__":

    def abbreviate(n: Optional[int], limit: int = 20) -> str:
        if n is None:
            return "None"

        s = str(n)
        if len(s) <= limit:
            return s

        half = (limit - 3) // 2
        return f"{s[:half]}...{s[-half:]}"

    def verify_case(a: int, b: int, c: int, x_expected: Optional[int],
                    y_expected: int):
        import time
        from datetime import timedelta

        parts = []
        if a != 0:
            parts.append(f"{abbreviate(a)}x^2")
        if b != 0:
            parts.append(f"{abbreviate(b)}x")
        parts.append(f"{abbreviate(c)}")
        rhs_str = " + ".join(parts).replace("+ -", "- ")

        print(f"Testing: y^2 = {rhs_str} ")
        start = time.time()
        result = KunerthSolver(a, b, c).solve()
        end = time.time()
        print(f"  ~  {timedelta(seconds=end-start)}")

        if result is None:
            if x_expected is None and y_expected == 0:
                print(f"  -> OK [Impossible]")
            else:
                print(f"  -> [FAIL] Exhausted search.")
            return

        x_calc, y_calc = result
        lhs = y_calc**2
        rhs = a * (x_calc**2) + b * x_calc + c

        if lhs != rhs:
            print(f"  -> [FAIL] Math Error!")
            return

        if a == 0:
            modular_root = abs(y_calc) % abs(b)
            status = "OK"
            note = "" if modular_root == y_expected else " (Found alternative root)"
            print(
                f"  -> {status}{note} [y: {abbreviate(y_calc)} , Mod: {abbreviate(modular_root)}]"
            )
        else:
            status = "OK" if x_calc == x_expected else f"OK (Found different solution)"
            print(
                f"  -> {status} [x: {abbreviate(x_calc)}, y: {abbreviate(y_calc)}]"
            )

    print("--- 1. Modular Cases (1878 Paper) ---")
    # Example 1: Page 341
    verify_case(0, 20521, 1124, None, 9758)
    # Example 2: Page 343
    verify_case(0, 21433, 16315, None, 10732)
    # Example 3: Page 343
    verify_case(0, 239, 127, None, 137)
    # Example 4: Page 343
    verify_case(0, 19139, 4163, None, 9179)

    print("--- 2. General Cases (1880 Paper) ---")
    # 1. Exempel: Page 358
    verify_case(61, 111, -101, -210, 1633)
    # 2. Exempel: Page 358
    verify_case(67, -290, 244, -24, 214)
    # 3. Exempel: Page 359
    verify_case(156, 447, 127, -41, 494)
    # 4. Exempel: Page 359
    verify_case(3245, 296, -5332, 67, 1767)
    # 5. Exempel: Page 360
    verify_case(639, -15, -615, -8, 201)
    # 6. Exempel: Page 361
    verify_case(126, 528, -264, -10, 84)
    # 7. Exempel: Page 361
    verify_case(540, -78, -110, -7, 164)
    # 8. Exempel: Page 362
    verify_case(199, 97, -20, -763, -10760)
    # 9. Exempel: Page 363
    verify_case(124, 37, 7, -787, 8762)
    # 10. Exempel: Page 363
    verify_case(46, -417, -324, 409, -2743)
    # 11. Exempel: Page 364
    verify_case(351, 358, -697, 43103, 807544)
    # 12. Exempel: Page 365
    verify_case(294, 371, 120, 8, 148)
    # 13. Exempel: Page 366
    verify_case(637, 0, 575, 1235, 31170)
    # 14. Exempel: Page 370
    verify_case(194, 63, 479, -83, 1154)
    # 15. Exempel: Page 371
    verify_case(290, 199, -123, -4, 61)

    print("--- 3. Impossible Cases ---")
    # Theoretical Basis: "Kriterien der Unmöglichkeit", 1880, p. 346
    # and Modular Exclusion, 1878b, p. 339.
    # y^2 = 8x + 5. Mod 8 check.
    verify_case(0, 8, 5, None, 0)
    # y^2 = 4x + 3. Mod 4 check.
    verify_case(0, 4, 3, None, 0)
    # y^2 = 100x + 2. Mod 4 check.
    verify_case(0, 100, 2, None, 0)