generate.py

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

3039 bytes, as of the version from 2026-09-23 03:43 (current). Recorded here, not run.

"""Lagrange numbers L_m of the Markov spectrum -- numberdb.org/T428

Run it with SageMath:

    $ sage -pip install numberdb          # once
    $ sage -python generate.py            # check the table against this code
    $ sage -python generate.py --publish  # send it, with NUMBERDB_API_KEY set

This table holds the values

    L_m = sqrt(9*m^2 - 4) / m

for Markov numbers m up to 10^12. The Markov numbers are obtained from the
Markov tree, starting at (1, 1, 1), using the Vieta involutions and normalising
each triple into increasing order.
"""

import os
import sys
from collections import deque

import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField


MARKOV_BOUND = 10**12
WORKING_GUARD = 64


def _key_from_stdin():
    if os.environ.get("NUMBERDB_KEY_FROM_STDIN") != "1":
        return
    token = sys.stdin.read().strip()
    if "=" in token and token.split("=", 1)[0].isupper():
        token = token.split("=", 1)[1].strip().strip("'\"")
    if token:
        os.environ["NUMBERDB_API_KEY"] = token


def normalised(triple):
    return tuple(sorted(ZZ(t) for t in triple))


def vieta_neighbours(triple):
    a, b, c = triple
    candidates = (
        (3*b*c - a, b, c),
        (a, 3*a*c - b, c),
        (a, b, 3*a*b - c),
    )
    for candidate in candidates:
        yield normalised(candidate)


def markov_triples(bound=MARKOV_BOUND):
    """Normalised Markov triples with largest entry at most bound."""
    root = (ZZ(1), ZZ(1), ZZ(1))
    seen = {root}
    pending = deque([root])
    while pending:
        triple = pending.popleft()
        yield triple
        for neighbour in vieta_neighbours(triple):
            if neighbour in seen or neighbour[-1] > bound:
                continue
            seen.add(neighbour)
            pending.append(neighbour)


def markov_numbers(bound=MARKOV_BOUND):
    return sorted({triple[-1] for triple in markov_triples(bound)})


def lagrange_number(m, digits):
    field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    m = ZZ(m)
    value = field(9*m*m - 4).sqrt() / field(m)
    if not value.is_finite():
        raise ArithmeticError("non-finite ball for m=%s" % (m,))
    return value


class LagrangeNumbersMarkovSpectrum(numberdb.Generator):
    table = os.environ.get("NUMBERDB_TABLE") or "T428"
    parameters = ("m",)
    type = "R"
    digits = 100
    rigour = "proven"

    def enumerate(self, bound=MARKOV_BOUND):
        for m in markov_numbers(bound):
            yield {"m": int(m)}

    def value(self, params, digits):
        return lagrange_number(params["m"], digits)


if __name__ == "__main__":
    _key_from_stdin()
    generator = LagrangeNumbersMarkovSpectrum()
    if os.environ.get("NUMBERDB_PUBLISH") == "1" or "--publish" in sys.argv:
        print(generator.publish(
            message="computed Lagrange numbers from Markov numbers"))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)