"""Laplacian energies of connected graphs -- numberdb.org/T548.

Run with SageMath:
    sage -pip install numberdb mpmath    # once
    sage -python generate.py             # verify every stored entry
    sage -python generate.py --fill      # fill the draft using NUMBERDB_API_KEY
    sage -python checks.py               # independent checks, no website writes

All connected simple graph isomorphism classes of orders 2 through 6 are
included, in order of vertex count and then ASCII graph6 identity. Identity
uses canonical_label(algorithm='sage'), validated with SageMath 10.9.
See argument-selection.md for the mathematical purpose and finite cutoff.

Integer Laplacian characteristic polynomials are factored over QQ. Their
roots in AA are compared exactly with the mean degree d=2m/n. The identity
LE=2*sum(mu-d for mu>d) avoids cancellation. Complete positive factors are
summed by their rational traces; partial factors use certified MPFI intervals.
The 80-bit guard gives 413 working bits at 100 digits, and precision doubles
if the width does not certify 105 relative decimal digits.

Exactness does not depend on observed digits. LE is rational precisely when
each irreducible factor lies entirely on one side of d. For the converse,
the chosen roots of each factor are its largest k roots. Every Galois
conjugate of their sum is at most that sum; when 0<k<degree, transitivity
moves the subset for some automorphism and makes it strictly smaller.
All factor multiplicities are positive, so contributions from other factors
cannot cancel that strict decrease. A rational total would be fixed by every
automorphism, a contradiction. This avoids constructing any compositum.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.graphs.graph import Graph
from sage.graphs.graph_generators import graphs
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.qqbar import AA
from sage.rings.real_mpfi import RealIntervalField

WORKING_GUARD = 80
MAX_ORDER = 6
EXPECTED_COUNTS = {2: 1, 3: 2, 4: 6, 5: 21, 6: 112}

@lru_cache(maxsize=None)
def graph_keys(max_order=MAX_ORDER):
    keys = []
    for n in range(2, max_order + 1):
        group = sorted({G.canonical_label(algorithm='sage').graph6_string()
                        for G in graphs(n) if G.is_connected()})
        if n in EXPECTED_COUNTS:
            assert len(group) == EXPECTED_COUNTS[n]
        keys.extend(group)
    return tuple(keys)


def graph_comment(G):
    n, m = G.order(), G.size()
    names = []
    if m == n * (n - 1) // 2:
        names.append(r'the complete graph $K_{%d}$' % n)
    degrees = G.degree()
    if m == n - 1 and max(degrees) <= 2:
        names.append(r'the path $P_{%d}$' % n)
    if n >= 3 and all(d == 2 for d in degrees):
        names.append(r'the cycle $C_{%d}$' % n)
    if G.is_bipartite():
        parts = G.bipartite_sets()
        a, b = sorted(map(len, parts))
        if m == a * b:
            names.append(r'the complete bipartite graph $K_{%d,%d}$' % (a,b))
    return ('This is ' + ', also '.join(names) + '.') if names else None


@lru_cache(maxsize=None)
def spectral_parts(key):
    G = Graph(key, format='graph6')
    assert G.order() >= 2 and G.is_connected()
    p = G.laplacian_matrix().charpoly().change_ring(QQ)
    d = QQ(2 * G.size()) / QQ(G.order())
    exact = QQ(0)
    partial = []
    for f, multiplicity in p.factor():
        roots = f.roots(ring=AA, multiplicities=False)
        assert len(roots) == f.degree()
        positive = [r for r in roots if r > d]
        if len(positive) == f.degree():
            trace = -f[f.degree()-1] / f.leading_coefficient()
            exact += 2 * multiplicity * (trace - f.degree() * d)
        elif positive:
            partial.extend((r, multiplicity) for r in positive)
    return exact, tuple(partial), d


def laplacian_energy(key, digits=100):
    exact, partial, d = spectral_parts(key)
    if not partial:
        return exact
    bits = numberdb.bits(digits, losing=WORKING_GUARD)
    while True:
        I = RealIntervalField(bits)
        result = I(exact)
        for root, multiplicity in partial:
            result += 2 * multiplicity * (I(root) - I(d))
        lower = result.lower().exact_rational()
        upper = result.upper().exact_rational()
        assert 0 < lower < upper
        if upper - lower < lower / ZZ(10)**(digits+5):
            return result
        bits *= 2
        if bits > 20000:
            raise ArithmeticError('Could not certify the requested precision')


class LaplacianEnergies(numberdb.Generator):
    table = 'T548'
    parameters = ('g',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'argument-selection.md',
             'census-reference.json', 'corpus-reference.json',
             'numerical-evidence.json', 'skill-evidence.json')

    def enumerate(self, max_order=MAX_ORDER):
        for key in graph_keys(max_order):
            yield {'g': key}

    def value(self, params, digits):
        key = str(params['g'])
        number = laplacian_energy(key, digits)
        comment = graph_comment(Graph(key, format='graph6'))
        return {'number': number, 'comment': comment} if comment else number


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=', 1)[0].strip().isidentifier():
            key = key.split('=', 1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = key
    generator = LaplacianEnergies()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill the complete 2 through 6 vertex census of Laplacian energies'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)