#!/usr/bin/env python3
"""Exact numerical certificates for arc_bound_q_two_thirds.pdf.

Run with Python 3.8 or later; there are no third-party dependencies:

    python3 verify_arc_bounds.py
    python3 verify_arc_bounds.py 128 512 2048

The output is one JSON object per field. All inequalities use integers.
This program searches necessary numerical conditions for an envelope
component. It does not construct or enumerate algebraic curves or arcs.

Section 8 of the accompanying write-up lists the conditions. The default
field sizes and their boundary certificates are asserted against the table.
"""

import argparse
import json
from math import isqrt


DEFAULT_FIELDS = (32, 128, 512, 2048, 8192, 32768, 131072, 524288)
# q: (unrefined first t, refined first t, nu, d, smallest D, largest D, max E)
EXPECTED = {
    32: (10, 10, 2, 10, 45, 45, 0),
    128: (24, 25, 4, 24, 136, 138, 8),
    512: (61, 62, 8, 61, 453, 457, 36),
    2048: (174, 175, 16, 174, 1876, 1881, 86),
    8192: (458, 459, 16, 458, 13079, 13081, 42),
    32768: (1009, 1010, 32, 1009, 31761, 31783, 720),
    131072: (2866, 2866, 64, 2866, 128208, 128298, 5778),
    524288: (7947, 7948, 64, 7947, 986632, 986669, 2414),
}


def ceildiv(a, b):
    """Ceiling of a/b, for integer a and positive integer b."""
    return -((-a) // b)


def restrict_quadratic(intervals, b, c):
    """Intersect closed integer intervals with z*z-b*z+c >= 0.

    Exact integer roots and nonsquare discriminants are both handled.
    Only the two outer intervals can survive.
    """
    discriminant = b * b - 4 * c
    if discriminant <= 0:
        return intervals
    left = (b - isqrt(discriminant)) // 2
    if left * left - b * left + c < 0:
        left -= 1
    right = b - left
    answer = []
    for start, stop in intervals:
        if start <= min(stop, left):
            answer.append((start, min(stop, left)))
        if max(start, right) <= stop:
            answer.append((max(start, right), stop))
    return answer


def passes_normal_form(q, order, degree, polar_base_degree):
    """The E<d criterion and its polynomial degree consequence."""
    if order < 4 or polar_base_degree >= degree:
        return True
    if (degree - 1) % order:
        return False
    if degree <= q and (order - 1) * degree > q - 1:
        return False
    return True


def admissible_parameters(q, t, normal_form=False):
    """Generate all tuples allowed by the necessary finite conditions.

    Call only in the range 3*t < q+5, t < q/2. The generator uses exact
    quadratic interval restrictions to avoid scanning irrelevant D values.
    """
    assert 3 * t < q + 5 and 2 * t < q
    k = q + 2 - t
    nu = 2
    while nu < t:
        mu = q // nu
        d_min = max(nu + 1, 3 + ceildiv(q + 2 - 2 * t, nu))
        for d in range(d_min, t + 1):
            # Local genus loss on C and the dual weighted point upper bound.
            D_min = max(
                mu + 1,
                ceildiv(d * (q + 2 * nu + 1 - d), nu * (nu + 1)),
                ceildiv((k + mu) * d, q + 1),
            )
            # Polar degree and N >= k*d on C.
            D_max = min(d * (d - 1) // nu, (t - 1) * d // nu)
            if D_min > D_max:
                continue
            intervals = [(D_min, D_max)]
            # D*(q+2*mu+1-D) <= mu*(mu+1)*d.
            intervals = restrict_quadratic(
                intervals, q + 2 * mu + 1, d * mu * (mu + 1)
            )
            # H(D) <= 0: dual smooth points versus original weighted count.
            intervals = restrict_quadratic(intervals, q + nu + 2, (q + 1) * d)
            # D*(D-1) >= mu*d.
            intervals = restrict_quadratic(intervals, 1, -mu * d)
            for start, stop in intervals:
                for D in range(start, stop + 1):
                    E = d * (d - 1) - nu * D
                    E_dual = D * (D - 1) - mu * d
                    if normal_form and not (
                        passes_normal_form(q, nu, d, E)
                        and passes_normal_form(q, mu, D, E_dual)
                    ):
                        continue
                    N_min = max(k * d, d * (q - d + 2), D * (q - D + 2))
                    N_max = min((q + 1) * d - nu * D, (q + 1) * D - mu * d)
                    if N_min <= N_max:
                        yield {
                            "nu": nu, "mu": mu, "d": d, "D": D,
                            "E": E, "E_dual": E_dual,
                            "N_min": N_min, "N_max": N_max,
                        }
        nu *= 2


def first_parameter(q, t, normal_form=False):
    return next(admissible_parameters(q, t, normal_form), None)


def threshold(q, normal_form=False):
    """First t allowed by the inequalities, using monotonicity in t."""
    lower, upper = 1, (q + 4) // 3
    assert first_parameter(q, lower, normal_form) is None
    if first_parameter(q, upper, normal_form) is None:
        raise ValueError("No admissible upper endpoint in the classical-exclusion range")
    while upper - lower > 1:
        middle = (lower + upper) // 2
        if first_parameter(q, middle, normal_form) is None:
            lower = middle
        else:
            upper = middle
    assert first_parameter(q, upper - 1, normal_form) is None
    example = first_parameter(q, upper, normal_form)
    assert example is not None
    return {"first_t": upper, "arc_upper_bound": q + 2 - upper,
            "admissible_numerical_tuple": example}


def scalar_threshold(q):
    """T(q) from Theorem 1.1, with an integer square root."""
    best = None
    minimizers = []
    nu = 2
    while nu < q:
        A = ceildiv(q + 3 * nu + 2, nu + 2)
        discriminant = (nu - 1) ** 2 + 4 * nu * (q + 2)
        B = (1 - nu + isqrt(discriminant)) // 2
        if B * (B + nu - 1) < nu * (q + 2):
            B += 1
        value = max(A, B)
        if best is None or value < best:
            best, minimizers = value, [nu]
        elif value == best:
            minimizers.append(nu)
        nu *= 2
    return {"T": best, "minimizing_orders": minimizers}


def validate_field(q):
    if q < 32 or q & (q - 1) or (q.bit_length() - 1) % 2 != 1:
        raise ValueError("q must be 2^h >= 32 with h odd")


def certify(q):
    validate_field(q)
    unrefined = threshold(q, False)
    refined = threshold(q, True)
    scalar = scalar_threshold(q)
    if q >= 512:
        assert unrefined["first_t"] >= scalar["T"]

    result = {"q": q, "duality": unrefined, "normal_form": refined,
              "scalar_optimization": scalar}
    if q in EXPECTED:
        t0, t1, nu, d, D_min, D_max, E_max = EXPECTED[q]
        assert unrefined["first_t"] == t0
        assert refined["first_t"] == t1
        boundary = list(admissible_parameters(q, t0))
        assert [(row["nu"], row["d"], row["D"]) for row in boundary] == [
            (nu, d, D) for D in range(D_min, D_max + 1)
        ]
        assert max(row["E"] for row in boundary) == E_max
        if q >= 512:
            assert unrefined["first_t"] == scalar["T"]
        if t1 > t0:
            assert nu >= 4 and E_max < d and (d - 1) % nu != 0
        result["complete_unrefined_boundary"] = {
            "nu": nu, "d": d, "D_min": D_min, "D_max": D_max,
            "E_max": E_max, "tuple_count": len(boundary),
        }
        result["writeup_table_verified"] = True
    return result


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("q", nargs="*", type=int, help="field sizes (default: eight fields in the write-up)")
    args = parser.parse_args()
    for q in args.q or DEFAULT_FIELDS:
        print(json.dumps(certify(q), sort_keys=True), flush=True)


if __name__ == "__main__":
    main()