"""Expected convex-hull vertex counts -- numberdb.org/T471.

Run with SageMath and the numberdb package:
    sage -pip install numberdb
    sage -python generate.py                  # verify all stored entries
    sage -python generate.py --fill           # fill the existing draft
    sage -python checks.py                    # independent local checks

Private draft reads and --fill need NUMBERDB_API_KEY. Alternatively set
NUMBERDB_KEY_FROM_STDIN=1 and pipe the key to stdin. Filling changes entries
only; this program has no publication or review operation.

The triangle and parallelogram expectations are exact rationals from Buchta
(1984), also displayed on slide 5 of his Rouen 2012 lecture. The disk formula
is Efron's expected-edge integral, transformed by t=-cos(theta). Since a
planar convex polygon has equally many edges and vertices, this counts vertices.
Uniformity is area measure on the filled body. Affine maps preserve the count.
The common range 4..20 contains the first nontrivial case and small samples.
The universal n=3 expectation is 3 and serves as a control, but is not stored.

An 80-bit guard is used with Arb's certified integration of an entire
trigonometric integrand. Finite, narrow enclosures are checked before return.
No midpoint is used to certify a value. The checks report the measured radii.
"""
import os
import sys
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField

GUARD_BITS = 80
REGIONS = ('equilateral-triangle', 'unit-square', 'unit-disk')


def triangle(n):
    return 2 * sum((QQ(1) / k for k in range(1, n)), QQ(0))


def parallelogram(n):
    s = sum(((1 - QQ(1) / 2**k) / k for k in range(1, n + 1)), QQ(0))
    return QQ(8) / 3 * (s - QQ(1) / (n * 2**n))


def ellipse(n, digits=100):
    C = ComplexBallField(numberdb.bits(digits, losing=GUARD_BITS))
    pi = C.pi()
    def integrand(u, analytic):
        theta = pi * u
        return theta.sin()**4 * (u - theta.sin()*theta.cos()/pi)**(n - 2)
    z = C(8*n*(n-1)) / 3 * C.integral(integrand, 0, 1)
    assert z.real().is_finite() and z.imag().is_finite()
    assert z.imag().contains_zero()
    value = z.real()
    assert value > 0 and value < n + 1
    R = value.parent()
    assert R(value.rad()) < abs(value) * R(10)**(-digits-10)
    return value


class ExpectedHullVertices(numberdb.Generator):
    table = 'T471'
    parameters = ('region', 'n')
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'dry_run.py')

    def enumerate(self):
        for region in REGIONS:
            for n in range(4, 21):
                yield {'region': region, 'n': n}

    def value(self, params, digits):
        n = int(params['n'])
        region = params['region']
        if region == 'equilateral-triangle':
            return triangle(n)
        if region == 'unit-square':
            return parallelogram(n)
        if region == 'unit-disk':
            ans = {'number': ellipse(n, digits)}
            if n == 4:
                ans['comment'] = r"The exact value is $4-35/(12\pi^2)$, obtained from the disk case of Sylvester's four-point problem."
            return ans
        raise ValueError('unsupported affine class')


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN'):
        numberdb.configure(api_key=sys.stdin.read().strip())
    g = ExpectedHullVertices()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_FILL'):
        print(g.publish(message='Fill 51 expected hull-vertex counts with exact rationals and certified disk integrals.', assisted_by='codex-cli'))
    else:
        report = g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)