back to table · edit · history · where entries came from · files · download
9805 bytes, as of the version from 2026-09-16 12:57 (current). Recorded here, not run.
"""Core thresholds of random $k$-uniform hypergraphs -- numberdb.org/T273.
For k >= 3 and r >= 2, c_{k,r} is the edge density m/n at which the random
k-uniform hypergraph with n vertices and m edges acquires a nonempty r-core.
Molloy and Cain-Wormald give
c_{k,r} = min over lambda > 0 of
lambda / (k * P(Poisson(lambda) >= r - 1)^(k - 1)).
This table stores c_{k,r} for 3 <= k <= 8 and 2 <= r <= 8 at 100 digits in
ball arithmetic.
Run it with SageMath:
$ sage -pip install numberdb # once
$ sage -python generate.py # check the table against this code
$ sage -python generate.py --publish # fill the draft, with NUMBERDB_API_KEY set
The minimiser lambda_{k,r} is enclosed by bisection on the sign of
g(lambda) = Q(lambda,r-1) - (k-1) exp(-lambda) lambda^(r-1)/(r-2)!,
where Q(lambda,s) = P(Poisson(lambda) >= s). The signs at both ends of the
final bracket are checked in ball arithmetic. Evaluating c_{k,r} on that
bracket gives the stored enclosure.
"""
import os
import sys
from decimal import Decimal
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfr import RealField
TABLE = os.environ.get("NUMBERDB_TABLE", "T273")
# Measured before filling the draft: 42 entries, longest value 119 characters,
# and an entries block of about 13 KB including comments.
K_MIN = 3
K_MAX = 8
R_MIN = 2
R_MAX = 8
# Bits of working precision beyond what the written digits need. Measured at
# 100 digits: every result ball has radius below 9e-106.
WORKING_GUARD = 64
# Decimal places of the bracket around the minimiser beyond the digits written.
# On this rectangle, the value enclosure is still under 9e-106 in radius.
BRACKET_GUARD = 6
COMMENT_DIGITS = 12
PUBLISHED_TEN_DIGITS = {
(3, 2): "0.8184691607",
(3, 3): "1.5528299396",
(4, 2): "0.7722798398",
(5, 2): "0.7017802665",
}
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 poisson_tail(x, s):
"""Q(x,s) = P(Poisson(x) >= s), for integer s >= 1 and a ball x."""
if s < 1:
raise ValueError("the finite tail formula here expects s >= 1")
total = x.parent()(0)
term = x.parent()(1)
for j in range(s):
if j:
term = term * x / j
total += term
return 1 - (-x).exp() * total
def minimiser_function(x, k, r):
"""The numerator of the derivative of lambda / Q(lambda,r-1)^(k-1)."""
return (
poisson_tail(x, r - 1)
- (k - 1) * (-x).exp() * x ** (r - 1) / ZZ(r - 2).factorial()
)
def threshold_from_lambda(x, k, r):
"""The hypergraph r-core threshold at a candidate minimiser."""
return x / (k * poisson_tail(x, r - 1) ** (k - 1))
def minimiser(k, r, digits):
"""A ball enclosing the unique positive minimiser lambda_{k,r}."""
if k < 2:
raise ValueError("the formula is used here only for k >= 2")
if r < 2:
raise ValueError("the formula is used here only for r >= 2")
bits = numberdb.bits(digits, losing=WORKING_GUARD)
RBF, RR = RealBallField(bits), RealField(bits)
delta = RR(10) ** (-(digits + BRACKET_GUARD))
lo = RR(r - 1) - RR(1) / RR(k - 1)
hi = RR(4 * (k + r) + 10)
if not (minimiser_function(RBF(lo), k, r) < 0):
raise ArithmeticError("left endpoint does not have negative sign for k=%s, r=%s" % (k, r))
while not (minimiser_function(RBF(hi), k, r) > 0):
hi *= 2
if hi > 10000:
raise ArithmeticError("right endpoint did not bracket the minimiser for k=%s, r=%s" % (k, r))
while hi - lo > delta / 4:
mid = (lo + hi) / 2
value = minimiser_function(RBF(mid), k, r)
if value.contains_zero():
raise ArithmeticError(
"could not decide minimiser sign at lambda=%s for k=%s, r=%s" % (mid, k, r)
)
if value < 0:
lo = mid
else:
hi = mid
x0 = (lo + hi) / 2
left = minimiser_function(RBF(x0 - delta), k, r)
right = minimiser_function(RBF(x0 + delta), k, r)
if not (left < 0 and right > 0):
raise ArithmeticError("no sign change across minimiser bracket for k=%s, r=%s" % (k, r))
return RBF(x0).add_error(delta)
def threshold(k, r, digits):
"""(c_{k,r}, lambda_{k,r}) as real balls."""
lam = minimiser(k, r, digits)
value = threshold_from_lambda(lam, k, r)
if not (0 < value and value < 100):
raise ArithmeticError("threshold came out as %s at lambda=%s for k=%s, r=%s" % (value, lam, k, r))
return value, lam
class HypergraphCoreThresholds(numberdb.Generator):
"""Generator for T273, the table of hypergraph core thresholds."""
table = TABLE
parameters = ("k", "r")
type = "R"
digits = 100
rigour = "proven"
def enumerate(self):
for k in range(K_MIN, K_MAX + 1):
for r in range(R_MIN, R_MAX + 1):
yield {"k": str(k), "r": str(r)}
def value(self, params, digits):
k = int(params["k"])
r = int(params["r"])
value, lam = threshold(k, r, digits)
comment = r"The minimum is attained at $\lambda_{%d,%d}=%s$." % (
k,
r,
to_text(lam, COMMENT_DIGITS),
)
return {"number": value, "comment": comment}
def _mpmath_threshold(k, r, digits=120):
"""Independent numerical value using mpmath's incomplete gamma function."""
from mpmath import mp
mp.dps = digits + 30
factorial = mp.factorial(r - 2)
def q(x):
return mp.gammainc(r - 1, 0, x, regularized=True)
def g(x):
return q(x) - (k - 1) * mp.e ** (-x) * x ** (r - 1) / factorial
left = mp.mpf(r - 1) - mp.mpf(1) / mp.mpf(k - 1)
right = mp.mpf(4 * (k + r) + 10)
while g(right) <= 0:
right *= 2
for _ in range(4 * digits):
mid = (left + right) / 2
if g(mid) < 0:
left = mid
else:
right = mid
root = (left + right) / 2
return root / (k * q(root) ** (k - 1))
def run_integrity_checks():
"""Checks that do not share the generator's finite-sum implementation."""
bits = numberdb.bits(120, losing=WORKING_GUARD)
RBF = RealBallField(bits)
tolerance = RBF("1e-95")
for k in range(K_MIN, K_MAX + 1):
for r in range(R_MIN, R_MAX + 1):
value, _ = threshold(k, r, 100)
independent_value = _mpmath_threshold(k, r)
independent = RBF(str(independent_value))
if not (abs(value - independent) < tolerance):
raise ArithmeticError("mpmath comparison failed for k=%d, r=%d" % (k, r))
expected = PUBLISHED_TEN_DIGITS.get((k, r))
if expected is not None:
places = len(expected.split(".")[1])
independent_text = str(independent_value)
rounded = Decimal(independent_text).quantize(Decimal("1e-%d" % places))
if expected is not None and not (independent_text.startswith(expected) or str(rounded) == expected):
raise ArithmeticError("published ten-digit check failed for k=%d, r=%d" % (k, r))
table = numberdb.table("T156")
graph_values = table["Numbers"]
for r in range(3, R_MAX + 1):
value, _ = threshold(2, r, 100)
stored_half = RBF(graph_values[str(r)]["number"]) / 2
if not (abs(value - stored_half) < tolerance):
raise ArithmeticError("T156 half-value comparison failed for r=%d" % r)
def fill_draft_once(generator, message):
"""Fill a fresh draft without the client's empty upsert probe."""
from numberdb._generate import (
_check_precision,
_check_rigour,
_producer,
_run_name,
_source_files,
)
from numberdb._write import Entries, attach, submit_entries, to_text
table = generator.table
run = _run_name(generator)
entries = Entries(*generator.parameters)
for params in generator.enumerate():
params = dict(params)
wanted = generator.digits_for(params)
entry = generator._entry(params, wanted)
value = entry["number"]
identity = ",".join(str(params[name]) for name in generator.parameters)
_check_rigour(generator, table, identity, value)
written = to_text(value, wanted, generator.format)
_check_precision(table, identity, written, wanted, lowering=False)
record = dict(entry)
record.pop("digits", None)
entries.add(**params, **record, digits=wanted)
answer = submit_entries(
table,
entries,
message=message,
produced_by=_producer(generator, os.environ.get("NUMBERDB_ASSISTED_BY", "")),
upsert=False,
run=run,
rigour=generator.rigour,
)
for name, body in sorted(_source_files(generator).items()):
attach(table, name, body, run=run, message=message, rigour=generator.rigour)
return answer
if __name__ == "__main__":
_key_from_stdin()
generator = HypergraphCoreThresholds()
run_integrity_checks()
if "--publish" in sys.argv or os.environ.get("NUMBERDB_PUBLISH") == "1":
print(fill_draft_once(
generator,
message=(
"random k-uniform hypergraph r-core thresholds c_{k,r} for "
"3 <= k <= 8 and 2 <= r <= 8"
),
))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)