"""CKSU two-factor wreath-product construction bounds -- numberdb.org/T662.

Run with SageMath:
    sage -pip install numberdb mpmath         # once
    sage -python generate.py                 # verify all stored entries
    sage -python checks.py                   # independent numerical checks
    sage -python generate.py --fill-draft    # write entries; requires API key

The last command fills the existing draft; it does not offer or publish it.
For m>=5, beta_m is the unique root in (2,3) of
    [2*m*(m-1)]**s = 2*m**3 + binomial(m**3,2)*2**s.
We bisect the strictly increasing equivalent equation
    [m*(m-1)]**s - binomial(m**3,2) - 2*m**3*2**(-s) = 0.
Each accepted endpoint has an Arb-certified sign. Exact rational endpoints
are retained, then enclosed in a ball. 96 extra working bits support
bisection to digits+24 bits; the retained accuracy is measured by checks.py.
The finite selection and the canonical-skill manifest are documented in the
attached argument-selection.md and numerical-evidence.md.
"""
import sys
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

WORKING_GUARD = 96
BRACKET_GUARD = 24
MIN_M, MAX_M = 5, 40


def certified_bisect(function, low, high, steps):
    low, high = QQ(low), QQ(high)
    assert function(low) < 0 and function(high) > 0
    for _ in range(steps):
        mid = (low + high) / ZZ(2)
        fmid = function(mid)
        assert fmid.is_finite()
        if fmid < 0:
            low = mid
        elif fmid > 0:
            high = mid
        else:
            raise ArithmeticError('Working precision did not certify a midpoint sign')
    assert function(low) < 0 and function(high) > 0
    return low, high


def beta_bracket(m, digits=100):
    m = ZZ(m)
    if m < MIN_M:
        raise ValueError('The construction-bound family is defined for integers m>=5')
    field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    q = m**3
    pairs = q*(q-1) // ZZ(2)
    log_u, log_two = field(m*(m-1)).log(), field(2).log()
    def residual(s):
        return (field(s)*log_u).exp() - pairs - 2*q*(-field(s)*log_two).exp()
    return certified_bisect(residual, 2, 3, numberdb.bits(digits, losing=BRACKET_GUARD))


def construction_bound(m, digits=100):
    low, high = beta_bracket(m, digits)
    field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    return field((low+high)/ZZ(2), (high-low)/ZZ(2))


class ConstructionBounds(numberdb.Generator):
    table = 'T662'
    parameters = ('m',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'argument-selection.md',
             'numerical-evidence.md', 'skill-manifest.json')

    def enumerate(self, max_m=MAX_M):
        for m in range(MIN_M, int(max_m)+1):
            yield {'m': str(m)}

    def value(self, params, digits):
        return construction_bound(ZZ(params['m']), digits)


if __name__ == '__main__':
    generator = ConstructionBounds()
    if '--fill-draft' in sys.argv:
        print(generator.publish(message='Fill 36 certified exact-spectrum construction bounds'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)