"""Specific heats of finite Ising rectangles at critical coupling -- numberdb.org/T593.

Run with SageMath:
    sage -pip install numberdb
    sage -python generate.py                 # verify every stored entry
    sage -python generate.py --fill          # fill existing draft with API key

The model is the simple vertex grid P_m square P_n, with open boundaries,
uniform ferromagnetic coupling, zero field, and K_c=log(1+sqrt(2))/2.
Values are C/(mn k_B), not unnormalized bond-energy variances.
All integer cut counts are recomputed; no attached numeric output is input.
See argument-selection.md for the discrete size family and finite extent.
The 96-bit guard is checked in checks.py, which records retained accuracy.
"""
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField
from functools import lru_cache
import os
import sys

MAX_SIDE = 6
GUARD_BITS = 96


@lru_cache(maxsize=None)
def cut_counts(m, n):
    """Coefficients of F(b)=sum_sigma b^(number of disagreeing bonds)."""
    m, n = int(m), int(n)
    if m < 1 or n < 1:
        raise ValueError('Side lengths must be positive integers')
    states = 1 << m
    inside = [sum(((s >> i) ^ (s >> (i + 1))) & 1
                  for i in range(m - 1)) for s in range(states)]
    current = [{inside[s]: 1} for s in range(states)]
    for _ in range(1, n):
        following = [{} for _ in range(states)]
        for s in range(states):
            target = following[s]
            for t in range(states):
                shift = inside[s] + (s ^ t).bit_count()
                for d, count in current[t].items():
                    target[d + shift] = target.get(d + shift, 0) + count
        current = following
    edges = 2*m*n-m-n
    counts = [0] * (edges + 1)
    for row in current:
        for d, count in row.items():
            counts[d] += count
    return tuple(ZZ(c) for c in counts)


def heat_from_counts(counts, m, n, digits):
    """Evaluate 4 K_c^2 Var(d)/(mn), with exact inputs and Arb enclosures."""
    rb = RealBallField(numberdb.bits(digits, losing=GUARD_BITS))
    b = rb(2).sqrt() - 1
    kc = (rb(1) + rb(2).sqrt()).log() / rb(2)
    moments = [rb(0), rb(0), rb(0)]
    power = rb(1)
    for d, count in enumerate(counts):
        weight = rb(count) * power
        moments[0] += weight
        moments[1] += d * weight
        moments[2] += (d*d) * weight
        power *= b
    variance = moments[2]/moments[0] - (moments[1]/moments[0])**2
    answer = rb(4)*kc**2*variance/rb(m*n)
    if not answer.is_finite() or not answer > 0:
        raise ArithmeticError('The heat must have a finite positive enclosure')
    return answer


class SpecificHeats(numberdb.Generator):
    table = 'T593'
    parameters = ('m', 'n')
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'argument-selection.md',
             'numerical-evidence.md', 'counts.json', 'checks-result.json',
             'skill-evidence.json')

    def enumerate(self, max_side=MAX_SIDE):
        for m in range(2, max_side + 1):
            for n in range(m, max_side + 1):
                yield {'m': m, 'n': n}

    def value(self, params, digits):
        m, n = int(params['m']), int(params['n'])
        if not 2 <= m <= n:
            raise ValueError('Require 2 <= m <= n')
        return heat_from_counts(cut_counts(m, n), m, n, digits)


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        token = sys.stdin.read().strip()
        if '=' in token and token.split('=', 1)[0].replace('_', '').isalnum():
            token = token.split('=', 1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = token
    g = SpecificHeats()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(g.publish(message='Fill 15 critical specific heats with certified 100-digit values',
                        assisted_by='codex-cli'))
    else:
        report = g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)