generate.py

back to table · edit · history · where entries came from · files · download

6609 bytes, as of the version from 2026-10-28 22:14. Recorded here, not run.

"""Lee–Yang zeros of finite critical Ising rectangles -- numberdb.org/T592.

Run with SageMath:
    sage -pip install numberdb              # once
    sage -python generate.py                 # verify every stored entry
    sage -python generate.py --fill-draft    # write entries, API key required

The write command fills T592; it never offers or publishes the table.
The default family is every 2 <= m <= n <= 6 with free boundaries. Exact
integer counts are computed by row transfer, then b = sqrt(2)-1 is substituted
in Q(sqrt(2)). Palindromy reduces the root problem to t = z + 1/z. All t roots
are isolated in AA and proved to lie in (-2,2); the two z roots are reconstructed
in certified balls. For odd mn the root -1 is separated exactly. Exact ordering
of real t roots, negative imaginary part first, is the family's lexicographic
ordering by (Re z, Im z). No rounded roots are used as input.

WORKING_GUARD is 112 bits beyond numberdb.bits(digits), checked by checks.py.
The count recurrence, root reconstruction, and exact -1 branch are described
in argument-selection.md and numerical-evidence.md. This source has no data
file dependencies; count-data.json is a reproducibility output only.
"""
import os
import sys
from functools import lru_cache

import numberdb.sage as numberdb
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.number_field.number_field import QuadraticField
from sage.rings.qqbar import AA, _init_qqbar
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField

# Named Sage imports do not run the qqbar initialisation normally done by
# Sage's aggregate import. Explicit initialisation supplies its exact constants.
_init_qqbar()
WORKING_GUARD = 112
SIDE_LIMIT = 6
F = QuadraticField(2, 'a', embedding=AA(2).sqrt())
BC = F.gen() - 1
R = PolynomialRing(F, 'z')
z = R.gen()


@lru_cache(maxsize=None)
def count_family(m, max_n=SIDE_LIMIT):
    """Return all integer coefficients c[k,d] for width m through max_n rows.

    A state is the last row's bit mask. Appending row s adds popcount(s)
    occupied sites and popcount(s xor old)+horizontal(s) cut bonds.
    """
    states = range(1 << m)
    pop = [s.bit_count() for s in states]
    horizontal = [((s ^ (s >> 1)) & ((1 << (m-1))-1)).bit_count() for s in states]
    dp = [{(pop[s], horizontal[s]): 1} for s in states]
    result = {}
    for n in range(1, max_n+1):
        if n >= m:
            total = {}
            for row in dp:
                for key, count in row.items():
                    total[key] = total.get(key, 0) + count
            result[n] = total
        if n == max_n:
            break
        nxt = []
        for s in states:
            row = {}
            for old in states:
                delta = pop[s ^ old] + horizontal[s]
                for (k, d), count in dp[old].items():
                    key = (k + pop[s], d + delta)
                    row[key] = row.get(key, 0) + count
            nxt.append(row)
        dp = nxt
    return result


def counts(m, n):
    return count_family(m, max(n, SIDE_LIMIT))[n]


@lru_cache(maxsize=None)
def field_polynomial(m, n):
    powers = [BC**d for d in range(2*m*n-m-n+1)]
    coeff = [F(0)] * (m*n+1)
    for (k, d), c in counts(m, n).items():
        coeff[k] += ZZ(c) * powers[d]
    return R(coeff)


@lru_cache(maxsize=None)
def root_data(m, n):
    """Exact real-part identities and multiplicities in the required order."""
    p = field_polynomial(m, n)
    assert p.degree() == m*n and p[0] == p[m*n] == 1
    assert p.list() == list(reversed(p.list()))
    # Extract -1 to handle the odd-degree self-reciprocal case exactly.
    minus_one = 0
    while p(-1) == 0:
        p, remainder = p.quo_rem(z+1)
        assert remainder == 0
        minus_one += 1
    assert p.degree() % 2 == 0
    assert p.list() == list(reversed(p.list()))
    h = p.degree() // 2
    T = PolynomialRing(F, 't'); t = T.gen()
    cheb = [T(2), t]
    for k in range(2, h+1):
        cheb.append(t*cheb[-1]-cheb[-2])
    q = T(p[h]) + sum((p[h+k]*cheb[k] for k in range(1,h+1)),T(0))
    # Independently check the reduction as a polynomial identity.
    reconstructed = R(0)
    for k,c in enumerate(q):
        reconstructed += c*z**(h-k)*(z*z+1)**k
    assert reconstructed == p
    real_roots = sorted(q.roots(ring=AA), key=lambda pair: pair[0])
    assert sum(int(e) for _, e in real_roots) == h
    assert all(-2 < t0 < 2 for t0, _ in real_roots)
    roots = []
    if minus_one:
        roots.append((AA(-1), AA(0), minus_one))
    for t0,e in real_roots:
        x = t0/ZZ(2)
        y = (AA(1)-x*x).sqrt()
        roots.extend([(x,-y,int(e)),(x,y,int(e))])
    assert sum(e for _,_,e in roots) == m*n
    return roots


def enclosed_root(m,n,j,digits):
    x,y,e = root_data(m,n)[j-1]
    if x == -1 and y == 0:
        return ZZ(-1),e
    rb = RealBallField(numberdb.bits(digits,losing=WORKING_GUARD))
    cb = ComplexBallField(rb.prec())
    value = cb(rb(x),rb(y))
    assert value.real().is_finite() and value.imag().is_finite()
    return value,e


class LeeYangZeros(numberdb.Generator):
    table = 'T592'
    parameters = ('m','n','j')
    type = 'C'
    digits = 100
    rigour = 'proven'

    def enumerate(self, limit=SIDE_LIMIT):
        for m in range(2,limit+1):
            for n in range(m,limit+1):
                for j in range(1,len(root_data(m,n))+1):
                    yield {'m':m,'n':n,'j':j}

    def value(self,params,digits):
        m,n,j = (int(params[k]) for k in self.parameters)
        value,mult = enclosed_root(m,n,j,digits)
        comment = 'Multiplicity %d.' % mult
        if value == -1:
            comment += ' The root is exactly $-1$; spin reversal makes every odd-degree field polynomial divisible by $z+1$.'
        return {'number':value,'comment':comment}


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].isidentifier():
            key = key.split('=',1)[1].strip().strip('"\'')
        numberdb.configure(api_key=key)
    generator = LeeYangZeros()
    if '--fill-draft' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill critical Ising Lee-Yang zeros from exact cut counts and certified algebraic roots',assisted_by='codex-cli'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)