"""Global minimum energies of the Thomson problem -- numberdb.org/T478.

Run with SageMath and the accompanying sources.json:
    sage -pip install numberdb
    sage -python generate.py
    sage -python generate.py --fill  # authenticated draft entries only

The default verifies every stored row. Set NUMBERDB_API_KEY to read a private
 draft; NUMBERDB_KEY_FROM_STDIN=1 alternatively reads the key from stdin.
Neither mode publishes or offers the draft for review.

E_n sums 1/Euclidean_distance over unordered pairs on the unit sphere.
Proven cases 2,3,4,5,6,12 use exact rational or radical expressions.
Other rows round Deeter's printed decimals to 12 significant digits, after
independent coordinate energy checks in checks.py. Normalising the rounded
coordinates exposes a last-digit discrepancy at n=22 and n=44 in the
13-digit source; twelve significant digits contain every recomputed value. Those checks do not prove global optimality.
The source snapshot and its dated manifest are attached. Only 2 <= n <= 100
is included, a consecutive range of small benchmark instances.
"""
from decimal import Decimal, localcontext
import json
import os
from pathlib import Path
import sys

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

WORKING_GUARD = 64  # checked against 200-digit independent distance sums
PROVEN = (2, 3, 4, 5, 6, 12)
SOURCES = json.loads(Path(__file__).with_name('sources.json').read_text())
COMMENTS = {
    2: r'Antipodal points give $E_2=1/2$. Optimality follows because a chord of the unit sphere has length at most $2$.',
    3: r'An equilateral triangle on a great circle gives $E_3=\sqrt3$. Optimality follows from the convexity bound in CITE{formula-simplex}.',
    4: r'The regular tetrahedron gives $E_4=3\sqrt6/2$. Optimality follows from the convexity bound in CITE{formula-simplex}.',
    5: r'The triangular bipyramid gives $E_5=1/2+3\sqrt2+\sqrt3$. Global optimality is proved in CITE{Schwartz}.',
    6: r'The regular octahedron gives $E_6=3/2+6\sqrt2$. Global optimality follows from CITE{CohnKumar}, Theorem 1.2.',
    12: r'The regular icosahedron gives $E_{12}=3+30/\sqrt{2-2/\sqrt5}+30/\sqrt{2+2/\sqrt5}$. Global optimality follows from CITE{CohnKumar}, Theorem 1.2.',
}


def closed_form(n, digits):
    if n == 2:
        return QQ(1) / 2
    R = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    if n == 3:
        return R(3).sqrt()
    if n == 4:
        return 3 * R(6).sqrt() / 2
    if n == 5:
        return R(QQ(1) / 2) + 3 * R(2).sqrt() + R(3).sqrt()
    if n == 6:
        return R(QQ(3) / 2) + 6 * R(2).sqrt()
    if n == 12:
        return 3 + 30 / (2 - 2 / R(5).sqrt()).sqrt() + 30 / (2 + 2 / R(5).sqrt()).sqrt()
    raise ValueError('No proven closed form implemented for n=%s' % n)


class ThomsonEnergies(numberdb.Generator):
    table = 'T478'
    parameters = ('n',)
    type = 'R'
    digits = 100
    rigour = 'heuristic (agreement-checked)'
    files = ('generate.py', 'checks.py', 'check.py', 'sources.json',
             'source-manifest.json', 'source-evidence.md', 'numerical-checks.json',
             'dry-run.txt', 'source-discrepancies.json',
             'stored-T478.json', 'stored-checks.json')

    def enumerate(self):
        for n in range(2, 101):
            yield {'n': str(n)}

    def value(self, params, digits):
        n = int(params['n'])
        if n in PROVEN:
            return {'number': closed_form(n, digits), 'comment': COMMENTS[n]}
        with localcontext() as context:
            context.prec = 12
            written = format(+Decimal(SOURCES['deeter'][str(n)]['energy']), 'f')
        return {'number': written, 'digits': 12}


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    generator = ThomsonEnergies()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        report = generator.publish(message='Fill the Thomson draft with independently checked source energies')
        print(report)
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)