"""Independent checks of Löbell growth rates -- numberdb.org/T612.

sage -python checks.py                # all local values and identities
sage -python checks.py --stored       # also the written decimal intervals

The optional key for a private draft is NUMBERDB_API_KEY. For the campaign,
NUMBERDB_KEY_FROM_STDIN=1 reads it from stdin; NUMBERDB_CHECK_STORED=1
selects the read-back check. No writes are made.

The dual graph is constructed from the two pentagonal flowers (Inoue §4),
then its cliques give finite special subgroups. Steinberg's formula constructs
the denominator without using the generator's radical formula. Its reciprocal
polynomial is bisected in exact rational arithmetic. Short element counts use
only involutions and commutations, independently of that rational series.
"""
import os
import sys
import json
import itertools
from decimal import Decimal
from fractions import Fraction
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.real_arb import RealBallField
from generate import LobellGrowth


def graph(n):
    # 0 and 1 are the bases. Upper faces: 2+i; lower: 2+n+i.
    adj = [set() for _ in range(2*n+2)]
    def edge(a, b):
        adj[a].add(b); adj[b].add(a)
    for i in range(n):
        u, v = 2+i, 2+n+i
        edge(0, u); edge(1, v)
        edge(u, 2+(i+1)%n); edge(v, 2+n+(i+1)%n)
        edge(u, v); edge(u, 2+n+(i-1)%n)
    return adj


def clique_counts(adj):
    counts = [1]
    cliques = [()]
    while cliques:
        extended = []
        for clique in cliques:
            start = clique[-1]+1 if clique else 0
            for v in range(start, len(adj)):
                if all(v in adj[u] for u in clique):
                    extended.append(clique+(v,))
        if extended:
            counts.append(len(extended))
        cliques = extended
    return counts


def canonical(word, adj):
    # Closure under all adjacent commutations and ss -> empty. This exhaustive
    # algorithm is only used on words of length at most 3.
    seen = {word}; pending = [word]
    while pending:
        w = pending.pop()
        for i in range(len(w)-1):
            if w[i] == w[i+1]:
                changed = w[:i]+w[i+2:]
            elif w[i+1] in adj[w[i]]:
                changed = w[:i]+(w[i+1],w[i])+w[i+2:]
            else:
                continue
            if changed not in seen:
                seen.add(changed); pending.append(changed)
    return min(seen, key=lambda w:(len(w),w))


def short_counts(adj):
    layer = {()}; counts = [1]
    for length in range(1,4):
        nxt = set()
        for w in layer:
            for s in range(len(adj)):
                reduced = canonical(w+(s,), adj)
                if len(reduced) == length:
                    nxt.add(reduced)
        counts.append(len(nxt)); layer = nxt
    return counts


def decimal_interval(text):
    d = Decimal(text)
    c = Fraction(d); e = d.as_tuple().exponent
    unit = Fraction(10**e) if e >= 0 else Fraction(1,10**(-e))
    return QQ(c.numerator)/c.denominator, QQ(unit.numerator)/unit.denominator


def bisect(poly, lo, hi, steps=440):
    lo,hi=QQ(lo),QQ(hi)
    assert poly(lo)<0<poly(hi)
    for _ in range(steps):
        middle=(lo+hi)/2
        if poly(middle)<0: lo=middle
        else: hi=middle
    assert poly(lo)<0<poly(hi)
    return lo,hi


def unit_power(n, R):
    rad=ZZ(n*(n-2)); d=rad.squarefree_part(); b=ZZ(rad//d).isqrt()
    assert d*b*b == rad and d>1
    # All units 1<epsilon<=tau have positive trace a <= trace(tau)+1;
    # a^2-d*b^2=+/-4, with the maximal-order parity conditions.
    candidates=[]
    for a in range(1,2*n+1):
        for c in range(1,2*n+1):
            if d%4 == 1:
                if (a-c)%2: continue
            elif a%2 or c%2: continue
            if abs(a*a-d*c*c) == 4:
                val=(R(a)+R(d).sqrt()*c)/2
                if val>1: candidates.append((val,QQ(a)/2,QQ(c)/2))
    assert candidates
    _,u,v=min(candidates,key=lambda row:row[0].mid())
    x,y=QQ(1),QQ(0)
    for power in range(1,100):
        x,y=x*u+d*y*v,x*v+y*u
        if (x,y)==(n-1,b):
            D=d if d%4==1 else 4*d
            return int(D),power,(u,v,int(d))
    raise AssertionError('No exact fundamental-unit power found')


def main(stored=False):
    # Controls before any use of the counting or bracketing procedures.
    assert short_counts([{1,2},{0,2},{0,1}]) == [1,3,3,1]
    assert short_counts([set(),set(),set()]) == [1,3,6,12]
    P=PolynomialRing(QQ,'t'); t=P.gen()
    lo,hi=bisect(t*t-2,1,2,80)
    assert lo*lo<2<hi*hi and hi-lo==QQ(1)/2**80
    print('Controls passed: (Z/2)^3, three free involutions, sqrt(2) bracket.',flush=True)
    R=RealBallField(500); g=LobellGrowth(); rows=[]; values={}
    doc=numberdb.table('T612') if stored else None
    regulators=numberdb.table('T131')['Numbers'] if stored else None
    if stored: assert list(doc['Numbers'])==[str(n) for n in range(5,31)]
    for params in g.enumerate():
        n=int(params['n']); result=g.value(params,100)
        ball=result['number'] if isinstance(result,dict) else result
        adj=graph(n); c=clique_counts(adj)
        assert c==[1,2*n+2,6*n,4*n]
        D=sum((-1)**k*c[k]*t**k*(1+t)**(3-k) for k in range(4))
        expected=(1-t)*(1-(2*n-2)*t+t*t)
        assert D==expected
        # x^3 D(1/x) is assembled by reversing coefficients.
        reciprocal=P(list(reversed(D.list())))
        lo,hi=bisect(reciprocal,2*n-3,2*n-2)
        assert ball.is_finite() and R(ball.rad())/R(ball.mid()) < R(10)**-117
        assert R(ball).upper()>=R(lo).lower() and R(ball).lower()<=R(hi).upper()
        q=ball*ball-(2*n-2)*ball+1
        assert q.is_finite() and q.contains_zero() and R(q.rad())<R(10)**-112
        assert (ball+1/ball-(2*n-2)).contains_zero()
        assert (ball.log()-R(n-1).arccosh()).contains_zero()
        assert ZZ(n*(n-2)).is_square() is False
        coeff=[]
        for k in range(4):
            coeff.append(int(((1+t)**3)[k])-sum(D[j]*coeff[k-j] for j in range(1,min(3,k)+1)))
        counts=short_counts(adj)
        assert coeff==counts
        disc,m,unit=unit_power(n,R)
        unit_ball=R(unit[0])+R(unit[1])*R(unit[2]).sqrt()
        assert (ball-unit_ball**m).contains_zero()
        text=to_text(ball,digits=100)
        if stored:
            entry=doc['Numbers'][str(n)]; text=entry['number']
            if str(disc) in regulators:
                rc,ru=decimal_interval(regulators[str(disc)]['number'])
                log_value=R(ball).log()/m
                assert R(rc-ru)<log_value<R(rc+ru)
        center,ulp=decimal_interval(text)
        left,right=center-ulp,center+ulp
        assert reciprocal(left)<0<reciprocal(right)
        # Containment, not merely overlap, for the serialized ball enclosure.
        assert R(left)<R(ball) < R(right)
        assert left<lo<hi<right
        values[n]=text
        rows.append({'n':n,'cliques':c,'short_counts':counts,'D':disc,'unit_power':m,
                     'relative_accuracy_bits':max(k for k in range(300, 430) if R(ball.rad())/R(ball.mid()) < R(2)**(-k)),'text':text})
        print('n=%d: incidence, short counts, rational bracket, unit power, serialization OK'%n,flush=True)
    print(json.dumps({'stored':stored,'rows':rows,'minimum_accuracy_bits':min(r['relative_accuracy_bits'] for r in rows)},indent=2))

if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    main('--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED')=='1')