"""Convex-position probabilities -- numberdb.org/T468.

Run with SageMath (or passagemath):
    sage -pip install numberdb
    sage -python generate.py                 # verify the stored draft/table
    sage -python generate.py --fill          # write entries with an authorized key

Set NUMBERDB_API_KEY for an authenticated draft read or fill. Filling entries
never publishes the draft. The campaign wrapper can instead supply the key on
stdin with NUMBERDB_KEY_FROM_STDIN=1 and select filling with NUMBERDB_PUBLISH=1.

The table covers n=3,...,8 for triangles, parallelograms and ellipses, using
uniform area measure, not boundary measure. The exact rational polygon values
use recurrences derived from Valtr's factorial formulas. The disk values use
Marckert, arXiv:1402.3512v1, equation (7), which tabulates 1-p, not p.
An ellipse has the disk's probability by invertible affine invariance.

An 80-bit guard beyond the requested precision is checked in checks.py, along
with exact factorial formulas, Dyck-path enumeration, a geometric area integral,
independent evaluation, and uniform-point simulations. The attached evidence
records the measured accuracy. No digits are inferred from a simulation.
"""
import os
import sys
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

GUARD_BITS = 80
REGIONS = ('equilateral-triangle', 'unit-square', 'unit-disk')
N_MIN, N_MAX = 3, 8
DISK_COMMENTS = {
    4: "This is Sylvester's four-point probability, $1-35/(12\\pi^2)$ CITE{Marckert}.",
    5: r'The exact value is $1-305/(48\pi^2)$ CITE{Marckert}.',
    6: r'The exact value is $1-(146400\pi^2-473473)/(11520\pi^4)$ CITE{Marckert}.',
    7: r'The exact value is $1-(512400\pi^2-2900611)/(23040\pi^4)$ CITE{Marckert}.',
    8: r'The exact value is $1-(62664108221+1721664000\pi^4-18670881600\pi^2)/(48384000\pi^6)$ CITE{Marckert}.',
}


def polygon_probability(region, n):
    """Exact recurrence from p(3)=1; all quotients are Sage rationals."""
    result = QQ(1)
    for j in range(3, n):
        if region == 'equilateral-triangle':
            result *= QQ(2*(3*j)*(3*j-1)*(3*j-2)) / QQ(j**3*(2*j+2)*(2*j+1))
        elif region == 'unit-square':
            result *= (QQ(2*(2*j-1)) / QQ(j*(j+1)))**2
        else:
            raise ValueError('not a polygon class in this table')
    return result


def disk_probability(n, digits):
    """Evaluate a polynomial in 1/pi^2 using ball arithmetic throughout."""
    if n == 3:
        return QQ(1)
    field = RealBallField(numberdb.bits(digits, losing=GUARD_BITS))
    q = field(1) / field.pi()**2
    if n == 4:
        result = 1 - field(QQ(35)/12)*q
    elif n == 5:
        result = 1 - field(QQ(305)/48)*q
    elif n == 6:
        result = 1 + q*(-field(QQ(146400)/11520) + q*field(QQ(473473)/11520))
    elif n == 7:
        result = 1 + q*(-field(QQ(512400)/23040) + q*field(QQ(2900611)/23040))
    elif n == 8:
        result = 1 + q*(-field(QQ(1721664000)/48384000)
                        + q*(field(QQ(18670881600)/48384000)
                             - q*field(QQ(62664108221)/48384000)))
    else:
        raise ValueError('the cited explicit disk formulas cover n=3,...,8')
    assert result.is_finite() and result > 0 and result < 1
    assert field(result.rad()) < abs(result) * field(10)**(-digits-10)
    return result


class ConvexPosition(numberdb.Generator):
    table = 'T468'
    parameters = ('region', 'n')
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'evidence-before-fill.json')

    def enumerate(self):
        for region in REGIONS:
            for n in range(N_MIN, N_MAX+1):
                yield {'region': region, 'n': n}

    def value(self, params, digits):
        n = int(params['n'])
        if not N_MIN <= n <= N_MAX:
            raise ValueError('n must be in the documented range 3,...,8')
        region = params['region']
        if region == 'unit-disk':
            result = {'number': disk_probability(n, digits)}
            if n in DISK_COMMENTS:
                result['comment'] = DISK_COMMENTS[n]
        elif region in REGIONS[:2]:
            result = {'number': polygon_probability(region, n)}
            if n == 4:
                result['comment'] = "This is Sylvester's four-point probability."
        else:
            raise ValueError('unknown region')
        if n == 3:
            result['comment'] = 'Three independent uniform area points are noncollinear with probability $1$.'
        return result


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    generator = ConvexPosition()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill 18 convex-position probabilities with exact rationals and certified disk enclosures'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)