"""Stationary points of the Hardy Z-function -- numberdb.org/T585.

Run with SageMath and the attached hardy_arb.pyx beside this file:
    sage -pip install numberdb
    sage -python generate.py                 # verify every stored entry
    sage -python checks.py                   # independent numerical checks

A working Sage Cython compiler and FLINT are required. To fill the authorized
private draft, pipe its key to stdin with NUMBERDB_KEY_FROM_STDIN=1 and
NUMBERDB_PUBLISH=1. This writes entries; it never offers or publishes a table.

The root index is certified by covering [0,B] with adjacent exact rational
intervals. On each interval either Z' excludes zero, or Z'' excludes zero
and the endpoint signs count exactly zero or one simple root. Subdivision
continues until one of these certificates succeeds. At the origin, evenness
and Z''>0 prove Z'(t)/t>0 on (0,1/1024]. No Riemann hypothesis is used.

After isolation, Newton iteration only proposes a location. Fresh strict
endpoint signs and Z'' exclusion prove its final enclosure at 493 bits for
100-digit output. The 160-bit guard and enclosure widths are measured in
numerical-evidence.json. See certification.md and argument-selection.md.
The bridge is reused from the attached source of numberdb.org/T584.
"""
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

class HardyZStationary(numberdb.Generator):
    table = 'T585'
    parameters = ('j',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py','hardy_arb.pyx','checks.py','check.py',
             'argument-selection.md','certification.md',
             'numerical-evidence.json','skill-manifest.json')

    def enumerate(self):
        for j in range(1,COUNT+1):
            yield {'j':str(j)}

    def value(self,params,digits):
        return root(int(params['j']),digits)

if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key=sys.stdin.read().strip()
        if '=' in key: key=key.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY']=key
    generator=HardyZStationary()
    if os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Store the certified first 100 positive Hardy Z stationary points'))
    else:
        report=generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)