"""Certify Hardy Z stationary arguments for numberdb.org/T586.

The isolation and refinement routines are reused from the live attached
T585 generator, fetched 2026-10-28. This module computes roots afresh;
no stored decimal argument is used. See certification.md for the proof.
"""
import os
import sys
import json
import hashlib
from functools import lru_cache
from pathlib import Path
import numberdb.sage as numberdb
from sage.misc.cython import cython_import
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ

COUNT = 100
COVER_BITS = 96
WORKING_GUARD = 160

@lru_cache(maxsize=1)
def bridge():
    return cython_import(str(Path(__file__).with_name('hardy_arb.pyx')))

def sign(x):
    if not x.is_finite():
        return 0
    return 1 if x > 0 else -1 if x < 0 else 0

def jet(a, b, prec, length=3):
    R = RealBallField(prec)
    x = R((QQ(a)+QQ(b))/2, (QQ(b)-QQ(a))/2)
    result = bridge().jet(x, length, prec)
    # The FLINT Hardy jet is real for real arguments. Nonfinite/wide balls
    # are never used as successful certificates; they cause subdivision.
    return [v.real() for v in result]

@lru_cache(maxsize=None)
def isolate(count=COUNT):
    roots = []
    leaves = []
    origin_end = QQ(1)/1024
    origin_second = 2*jet(0, origin_end, COVER_BITS)[2]
    assert sign(origin_second) == 1
    leaves.append(('0', str(origin_end), 'origin', 1))

    def cover(a, b, depth=0):
        coefficients = jet(a, b, COVER_BITS)
        first_sign = sign(coefficients[1])
        if first_sign:
            leaves.append((str(a), str(b), 'exclude', first_sign))
            return
        second_sign = sign(coefficients[2])
        if second_sign:
            sa = sign(jet(a, a, COVER_BITS, 2)[1])
            sb = sign(jet(b, b, COVER_BITS, 2)[1])
            if sa and sb:
                if sa != sb:
                    assert sb == second_sign and sa == -second_sign
                    roots.append((a,b))
                    leaves.append((str(a), str(b), 'root', second_sign))
                else:
                    leaves.append((str(a), str(b), 'monotone-exclude', sa))
                return
        if depth >= 30:
            raise ArithmeticError('Unresolved interval: %s, %s' % (a,b))
        midpoint = (a+b)/2
        cover(a,midpoint,depth+1)
        cover(midpoint,b,depth+1)

    a = origin_end
    k = 0
    while len(roots) < count:
        k += 1
        if k > 4000:
            raise ArithmeticError('Increase the explicit height budget')
        b = QQ(k)/4
        cover(a,b)
        a = b
    assert len(roots) == count
    # Check exact adjacency, not just the total length of the cover.
    cursor = QQ(0)
    totals = {}
    for left,right,method,s in leaves:
        assert QQ(left) == cursor and QQ(right) > QQ(left)
        cursor = QQ(right)
        totals[method] = totals.get(method,0)+1
    assert cursor == a and totals['root'] == count
    encoded = json.dumps(leaves,separators=(',',':')).encode()
    summary = {'count':count,'endpoint':str(a),'cover_bits':COVER_BITS,
               'leaf_counts':totals,'leaves_sha256':hashlib.sha256(encoded).hexdigest(),
               'origin_second_derivative':str(origin_second),
               'origin_endpoint':str(origin_end)}
    return tuple(roots), summary, leaves

@lru_cache(maxsize=None)
def root(j, digits=100):
    if not 1 <= j <= COUNT:
        raise ValueError('This certified implementation covers j=1,...,100')
    (a,b) = isolate()[0][j-1]
    prec = numberdb.bits(digits,losing=WORKING_GUARD)
    R = RealBallField(prec)
    x = R((a+b)/2)
    # All division here is in the real ball field. Midpoints are used only
    # as candidates, never as a claim of exactness for the root.
    for _ in range(18):
        values = bridge().jet(x,3,prec)
        f = values[1].real()
        df = 2*values[2].real()
        assert sign(df)
        new = x-f/df
        assert new.is_finite()
        x = R(new.mid())
        assert R(a) < x and x < R(b)
    centre = QQ(x.mid())
    radius = QQ(1) / (2**numberdb.bits(digits+15,losing=0))
    left, right = centre-radius, centre+radius
    assert a < left < right < b
    fl = jet(left,left,prec,2)[1]
    fr = jet(right,right,prec,2)[1]
    second = 2*jet(left,right,prec)[2]
    assert sign(fl)*sign(fr) == -1
    assert sign(second) == sign(fr)
    result = R(centre,radius)
    assert result.is_finite() and result.accuracy() >= numberdb.bits(digits)+24
    return result