"""Best known packings of equal circles in a square -- numberdb.org/T477.

Run with SageMath, after downloading this file and square-packings-source.json
from the table's attached files into the same directory:

    sage -pip install numberdb
    sage -python generate.py                  # compare every stored entry
    sage -python checks.py                    # independent coordinate checks
    sage -python generate.py --fill           # write entries, never publish

The last command needs NUMBERDB_API_KEY for a writable draft. The build runner
may instead pipe the key with NUMBERDB_KEY_FROM_STDIN=1. No action here changes
publication or review state.

The source is Packomania's 10-Sep-2026 index and its coordinate files retrieved
on 2026-10-25. Each JSON row preserves the printed radius, distance, side length
and all coordinates, together with the coordinate file's SHA-256 hash.
We use the HTML index, since the downloadable distance.txt has a different
last few digits in the two-circle row. General rows retain 10 significant
digits, selected by checks.py from agreement with independently evaluated
coordinates and all three source columns. These are approximations to known
configurations, not enclosures of an unknown globally optimal value.

Exact radical expressions use Arb at 100 digits plus 64 guard bits; checks.py
measures the resulting widths and repeats at 160 digits. The five square-grid
cases are exact rationals. The certified separation intervals for n=31,32,33
are from Markot (2021), section 7.5, DOI 10.1007/s10898-021-01086-z; their
converted values are written to 14 significant digits. All other closed forms
are given in the entry comments with their source.
"""
import json
import os
import sys
from decimal import Decimal, localcontext
from fractions import Fraction
from pathlib import Path

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

GUARD_BITS = 64
SNAPSHOT = json.loads(Path(__file__).with_name('square-packings-source.json').read_text())
ROWS = SNAPSHOT['rows']
GRID = {4: 2, 9: 3, 16: 4, 25: 5, 36: 6}
FORMULAS = {
    2: r'2+\sqrt2', 3: r'2+(\sqrt2+\sqrt6)/2',
    5: r'2+2\sqrt2', 6: r'2+12/\sqrt{13}', 7: r'4+\sqrt3',
    8: r'2+\sqrt2+\sqrt6', 12: r'2+15\sqrt{2/17}',
    14: r'6+\sqrt3', 15: r'4+\sqrt2+\sqrt6',
    18: r'2+24/\sqrt{13}', 20: r'(130+16\sqrt2)/17',
    23: r'2+2\sqrt2+2\sqrt6', 24: r'6+\sqrt2+\sqrt6',
    27: r'2+80/\sqrt{89}', 30: r'(252+10\sqrt{10})/26',
    39: r'2+30\sqrt{2/17}', 42: r'(434+24\sqrt3)/37',
    52: r'2+168/\sqrt{193}', 56: r'(688+14\sqrt{14})/50',
    99: r'2+340/\sqrt{389}',
}
CERTIFIED = {
    31: ('0.2175472916191243', '0.2175472916191280'),
    32: ('0.2131745625898764', '0.2131745625898779'),
    33: ('0.2113283841432631', '0.2113283841432645'),
}


def rational(text):
    """Parse decimal digits exactly, without QQ's optional decimal parser."""
    f = Fraction(str(text))
    return QQ(f.numerator) / QQ(f.denominator)


def side_exact(n, digits=100):
    if n in GRID:
        return QQ(2 * GRID[n])
    b = RealBallField(numberdb.bits(digits, losing=GUARD_BITS))
    root = lambda k: b(k).sqrt()
    if n == 2: return b(2) + root(2)
    if n == 3: return b(2) + (root(2) + root(6)) / 2
    if n == 5: return b(2) + 2 * root(2)
    if n == 6: return b(2) + b(12) / root(13)
    if n == 7: return b(4) + root(3)
    if n == 8: return b(2) + root(2) + root(6)
    if n == 12: return b(2) + 15 * (b(2) / 17).sqrt()
    if n == 14: return b(6) + root(3)
    if n == 15: return b(4) + root(2) + root(6)
    if n == 18: return b(2) + b(24) / root(13)
    if n == 20: return (b(130) + 16 * root(2)) / 17
    if n == 23: return b(2) + 2 * root(2) + 2 * root(6)
    if n == 24: return b(6) + root(2) + root(6)
    if n == 27: return b(2) + b(80) / root(89)
    if n == 30: return (b(252) + 10 * root(10)) / 26
    if n == 39: return b(2) + 30 * (b(2) / 17).sqrt()
    if n == 42: return (b(434) + 24 * root(3)) / 37
    if n == 52: return b(2) + b(168) / root(193)
    if n == 56: return (b(688) + 14 * root(14)) / 50
    if n == 99: return b(2) + b(340) / root(389)
    raise ValueError('No exact formula declared for n=%d' % n)


def rounded(text, digits=10):
    with localcontext() as ctx:
        ctx.prec = digits
        return format(+Decimal(text), 'f')


def proof_comment(n):
    if n <= 20:
        return 'Optimality is proven; see CITE{Peikert} and CITE{Nurmela}.'
    if n <= 27:
        return 'Optimality is proven in CITE{Nurmela}.'
    if n <= 30:
        return 'Optimality is proven by the interval computations in CITE{MarkotCsendes} and the structural verification in CITE{Markot2007}.'
    if n in CERTIFIED:
        return 'The global optimum is enclosed by the certified interval in CITE{Markot2021}, Section 7.5.'
    if n == 36:
        return 'Optimality is proven in CITE{KirchnerWengerodt}.'
    return 'This is a best known configuration; global optimality is not asserted.' if n in FORMULAS else ''


class SquarePackings(numberdb.Generator):
    table = 'T477'
    parameters = ('n', 'normalisation')
    type = 'R'
    digits = 100
    rigour = 'heuristic (agreement-checked)'
    files = ('generate.py', 'checks.py', 'square-packings-source.json')

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

    def value(self, params, digits):
        n = int(params['n'])
        norm = params['normalisation']
        if norm not in ('separation', 'container'):
            raise ValueError(norm)
        comment = proof_comment(n)
        if n in GRID or n in FORMULAS:
            side = side_exact(n, digits)
            number = side if norm == 'container' else QQ(2) / (side - 2)
            if n in GRID:
                k = GRID[n]
                fact = r'The $%d\times%d$ square grid gives $S_{%d}=%d$ and $d_{%d}=1/%d$.' % (k, k, n, 2*k, n, k-1)
            else:
                source = 'Friedman' if n <= 24 else 'Cantrell'
                fact = r'The exact container side is $S_{%d}=%s$ CITE{%s}; $d_{%d}=2/(S_{%d}-2)$.' % (n, FORMULAS[n], source, n, n)
            return {'number': number, 'comment': ' '.join(x for x in (fact, comment) if x)}
        if n in CERTIFIED:
            field = RealIntervalField(numberdb.bits(100, losing=GUARD_BITS))
            lo, hi = CERTIFIED[n]
            d = field(rational(lo), rational(hi))
            number = d if norm == 'separation' else 2 + 2 / d
            return {'number': number, 'digits': 14, 'comment': comment}
        number = rounded(ROWS[str(n)][norm], 10)
        record = {'number': number, 'digits': 10}
        if comment:
            record['comment'] = comment
        return record


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    generator = SquarePackings()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill checked square packing snapshot', assisted_by=os.environ.get('NUMBERDB_ASSISTED_BY', '')))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)