back to table · edit · history · where entries came from · files · download
17903 bytes, as of the version from 2026-09-09 09:14 (current). Recorded here, not run.
"""Periodic windows of the logistic map.
For f_r(x) = r*x*(1 - x), this generator stores, by superstable kneading word,
the window onset, superstable parameter, first period-doubling point, and
topological entropy for the primitive windows through period 8, excluding the
main period-doubling cascade.
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
"""
import os
import sys
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
DIGITS = 100
# Measured at 100 digits: recomputing all endpoint and entropy rows at 760 and
# 980 guard bits gives the same 100 written digits. The defining-equation
# residuals are below 1e-180 away from the two harmonic onsets, where the value
# is copied from the parent's period-doubling endpoint because the multiplier
# +1 continuation is singular.
WORKING_GUARD = 760
WORKING_BITS = numberdb.bits(DIGITS, losing=WORKING_GUARD)
RR = RealField(WORKING_BITS)
WINDOW_DATA = (
("RLC", 3, "3.83187405528331556841036277549610655579782785260369463047889"),
("RLLC", 4, "3.96027012722115260157125779850142054454978143186380426963011"),
("RLRRC", 5, "3.73891491297068499308521289037628945123011655080787136700626"),
("RLLRC", 5, "3.90570646983129035482040709205270692520123591507403554654039"),
("RLLLC", 5, "3.99026704697370151827752638855431798336598982156607813701226"),
("RLRRRC", 6, "3.62755752951552320337487410201744167984471413668637428601269"),
("RLLRLC", 6, "3.84456879219443297290078487624558310353185746845492078015162"),
("RLLRRC", 6, "3.93753644475455110456195058197103873365215915515224038284069"),
("RLLLRC", 6, "3.97776642226547291402054545126087675864272268564042588450143"),
("RLLLLC", 6, "3.99758311825456726610019472081786309367501729922294434876419"),
("RLRRRRC", 7, "3.70176915353795611419331084301300604849413321574933283184849"),
("RLRRLRC", 7, "3.77421418890091323212739972204104818206084012028500672645025"),
("RLLRLRC", 7, "3.88604587818782201152249198516686183130618919773650386210825"),
("RLLRRRC", 7, "3.92219340330970002847241358035937377280274034284494160144631"),
("RLLRRLC", 7, "3.95103216476130619789149231869205357950097805625282226150577"),
("RLLLRLC", 7, "3.96897685695553797149529480553910894603778370379935112154288"),
("RLLLRRC", 7, "3.98474761881553890777227690756145467967695741149867955409540"),
("RLLLLRC", 7, "3.99453780911119719041722057353600587626551227740194151958958"),
("RLLLLLC", 7, "3.99939706096209841119149570748595344287447837974184760843867"),
("RLRRRRRC", 8, "3.66219250368657675340876010432325910893615015663642213835629"),
("RLRRLRRC", 8, "3.80077094387466976937368109454104253124347270335887952260225"),
("RLLRLRRC", 8, "3.87054098436375722318755269464027944382282435049289399007371"),
("RLLRLRLC", 8, "3.89946895097060233416934344787695205856357586340214333515350"),
("RLLRRRLC", 8, "3.91204662107512651104763023851459235523482847662900300073974"),
("RLLRRRRC", 8, "3.93047299573214484064402135589651638257113463053635152580303"),
("RLLRRLRC", 8, "3.94421349573255098808610901732878288595325948829461131171778"),
("RLLLRLLC", 8, "3.96093369754758146637114229913583917292042387207979678535427"),
("RLLLRLRC", 8, "3.97372425567498033782288399372111449028970594265243660107476"),
("RLLLRRRC", 8, "3.98140895441526149539552930197501172904253702383791585384259"),
("RLLLRRLC", 8, "3.98774549537007548382140455558237811663416743658051704703934"),
("RLLLLRLC", 8, "3.99251952328433189778669097627568608641243470578643252758917"),
("RLLLLRRC", 8, "3.99621959590116485491127851415714485231231646692996928535121"),
("RLLLLLRC", 8, "3.99864163620644360240349251525825508554511035102421802984503"),
("RLLLLLLC", 8, "3.99984936201385107048028702522067059113459005561707696664198"),
)
WINDOWS = tuple(row[0] for row in WINDOW_DATA)
PERIOD = {word: period for word, period, _seed in WINDOW_DATA}
SUPERSTABLE_SEEDS = {word: seed for word, _period, seed in WINDOW_DATA}
#: What the table holds for each window, in the order it is shown.
#:
#: The three parameter values are held in both normalisations, because a
#: reader who works in $z^2+c$ should not have to convert: c was in the entry
#: comments, which is not a place a value can be searched for or cited.
#:
#: One parameter with seven values rather than a `normalisation` parameter
#: beside a `quantity` one, as T168 has: T168's grid is a rectangle and this
#: is not. The topological entropy has no normalisation -- it is a property of
#: the map, the same number whichever coordinate names it -- and a
#: normalisation parameter would have to claim either that it has both or
#: that it has one.
QUANTITIES = (
"onset", "onset-c",
"superstable", "superstable-c",
"doubling", "doubling-c",
"entropy",
)
#: The `r` row each `c` row converts.
BASE_OF = {q: q[:-2] for q in QUANTITIES if q.endswith("-c")}
#: The one row a theorem gives exactly. The period-3 window opens at
#: r = 1 + 2 sqrt 2, so c = -r(r-2)/4 = -(1 + 2 sqrt 2)(2 sqrt 2 - 1)/4
#: = -(8 - 1)/4 and the surd cancels. A decimal would say a number known to be
#: -7/4 is known to a hundred places. Checked against the solve below.
EXACT = {
("RLC", "onset-c"): "-7/4",
}
HARMONIC_ONSETS = {
"RLLRLC": "RLC",
"RLLLRLLC": "RLLC",
}
ENTROPY_POLYNOMIALS = {
"RLC": "x^2 - x - 1",
"RLLC": "x^3 - x^2 - x - 1",
"RLRRC": "x^4 - x^3 - x^2 + x - 1",
"RLLRC": "x^4 - x^3 - x^2 - x + 1",
"RLLLC": "x^4 - x^3 - x^2 - x - 1",
"RLRRRC": "x^4 - x^2 - 1",
"RLLRLC": "x^2 - x - 1",
"RLLRRC": "x^5 - x^4 - x^3 - x^2 + x - 1",
"RLLLRC": "x^4 - 2*x^3 + x^2 - 2*x + 1",
"RLLLLC": "x^5 - x^4 - x^3 - x^2 - x - 1",
"RLRRRRC": "x^3 - x^2 - 1",
"RLRRLRC": "x^6 - x^5 - x^4 + x^3 - x^2 - x + 1",
"RLLRLRC": "x^6 - x^5 - x^4 - x^3 + x^2 + x - 1",
"RLLRRRC": "x^3 - 2*x^2 + x - 1",
"RLLRRLC": "x^6 - x^5 - x^4 - x^3 + x^2 - x - 1",
"RLLLRLC": "x^6 - x^5 - x^4 - x^3 - x^2 + x + 1",
"RLLLRRC": "x^6 - x^5 - x^4 - x^3 - x^2 + x - 1",
"RLLLLRC": "x^6 - x^5 - x^4 - x^3 - x^2 - x + 1",
"RLLLLLC": "x^6 - x^5 - x^4 - x^3 - x^2 - x - 1",
"RLRRRRRC": "x^6 - x^4 - x^2 - 1",
"RLRRLRRC": "x^7 - x^6 - x^5 + x^4 - x^3 - x^2 + x - 1",
"RLLRLRRC": "x^5 - x^4 - 2*x^2 + x - 1",
"RLLRLRLC": "x^7 - x^6 - x^5 - x^4 + x^3 + x^2 - x - 1",
"RLLRRRLC": "x^6 - x^4 - 2*x^3 - x^2 - 2*x - 1",
"RLLRRRRC": "x^7 - x^6 - x^5 - x^4 + x^3 - x^2 + x - 1",
"RLLRRLRC": "x^7 - x^6 - x^5 - x^4 + x^3 - x^2 - x + 1",
"RLLLRLLC": "x^3 - x^2 - x - 1",
"RLLLRLRC": "x^7 - x^6 - x^5 - x^4 - x^3 + x^2 + x - 1",
"RLLLRRRC": "x^7 - x^6 - x^5 - x^4 - x^3 + x^2 - x + 1",
"RLLLRRLC": "x^6 - 2*x^5 + x^4 - 2*x^3 + x^2 - 1",
"RLLLLRLC": "x^7 - x^6 - x^5 - x^4 - x^3 - x^2 + x + 1",
"RLLLLRRC": "x^5 - x^4 - 2*x^3 + x - 1",
"RLLLLLRC": "x^6 - 2*x^5 + x^4 - 2*x^3 + x^2 - 2*x + 1",
"RLLLLLLC": "x^7 - x^6 - x^5 - x^4 - x^3 - x^2 - x - 1",
}
SPECIAL_ENDPOINT_COMMENTS = {
("RLC", "onset"): r"$r_{\mathrm{on}}=1+2\sqrt2$ CITE{OEISA086178}",
("RLC", "superstable"): (
r"$r_W$ is the root near $3.83187$ of "
r"$r^6-6r^5+4r^4+24r^3-16r^2-32r-64$"
),
("RLC", "doubling"): (
r"$r_{\mathrm{pd}}$ is the root near $3.84150$ of "
r"$r^6-6r^5+4r^4+24r^3-14r^2-36r-81$ CITE{OEISA086179}"
),
("RLRRC", "onset"): (
r"$r_{\mathrm{on}}$ is the first period-$5$ onset listed in "
r"OEIS A118452 CITE{OEISA118452}"
),
("RLRRRC", "onset"): (
r"$r_{\mathrm{on}}$ is the first period-$6$ onset listed in "
r"OEIS A118453 CITE{OEISA118453}"
),
("RLRRRRC", "onset"): (
r"$r_{\mathrm{on}}$ is the first period-$7$ onset listed in "
r"OEIS A118746 CITE{OEISA118746}"
),
}
_CACHE = {}
def set_working_bits(bits):
"""Use another RealField for agreement checks."""
global WORKING_GUARD, WORKING_BITS, RR, _CACHE
WORKING_GUARD = max(0, int(bits) - numberdb.bits(DIGITS, losing=0))
WORKING_BITS = int(bits)
RR = RealField(WORKING_BITS)
_CACHE = {}
def decimal(value, digits=DIGITS):
return RR(value).str(digits=int(digits))
def f_value(r, x):
return r * x * (1 - x)
def c_parameter(r):
r = RR(r)
return -r * (r - 2) / 4
def critical_orbit_with_dr(r, steps):
y = RR(1) / RR(2)
dr = RR(0)
for _ in range(steps):
dr = y * (1 - y) + r * (1 - 2 * y) * dr
y = f_value(r, y)
return y, dr
def orbit_derivatives(x, r, steps):
"""f_r^steps(x) and derivatives for the multiplier equation."""
y = x
dx = RR(1)
dr = RR(0)
dxx = RR(0)
dxr = RR(0)
for _ in range(steps):
a = 1 - 2 * y
y_next = r * y * (1 - y)
dx_next = r * a * dx
dr_next = y * (1 - y) + r * a * dr
dxx_next = r * (a * dxx - 2 * dx * dx)
dxr_next = a * dx + r * (a * dxr - 2 * dr * dx)
y, dx, dr, dxx, dxr = y_next, dx_next, dr_next, dxx_next, dxr_next
return y, dx, dr, dxx, dxr
def superstable_parameter(word, digits=DIGITS):
key = ("superstable", word, WORKING_BITS)
if key in _CACHE:
return _CACHE[key]
p = PERIOD[word]
r = RR(SUPERSTABLE_SEEDS[word])
target = RR(10) ** (-(int(digits) + 90))
for _ in range(80):
y, dr = critical_orbit_with_dr(r, p)
step = (y - RR(1) / RR(2)) / dr
r -= step
if abs(step) < target:
break
y, _ = critical_orbit_with_dr(r, p)
if abs(y - RR(1) / RR(2)) > target:
raise ArithmeticError("%s: superstable equation not solved" % word)
_CACHE[key] = r
return r
def solve_multiplier(word, target_mu, digits=DIGITS):
key = ("multiplier", word, int(target_mu), WORKING_BITS)
if key in _CACHE:
return _CACHE[key]
p = PERIOD[word]
r = superstable_parameter(word, digits)
x = RR(1) / RR(2)
target = RR(10) ** (-(int(digits) + 90))
for j in range(1, 81):
mu = RR(target_mu) * RR(j) / RR(80)
for _ in range(40):
y, dx, dr, dxx, dxr = orbit_derivatives(x, r, p)
f1 = y - x
f2 = dx - mu
j11 = dx - 1
j12 = dr
j21 = dxx
j22 = dxr
det = j11 * j22 - j12 * j21
step_x = (f1 * j22 - j12 * f2) / det
step_r = (j11 * f2 - f1 * j21) / det
x -= step_x
r -= step_r
if max(abs(step_x), abs(step_r)) < target:
break
y, dx, _, _, _ = orbit_derivatives(x, r, p)
if abs(y - x) > target or abs(dx - RR(target_mu)) > target:
raise ArithmeticError("%s: multiplier %s equation not solved" % (
word, target_mu))
_CACHE[key] = r
return r
def onset_parameter(word, digits=DIGITS):
parent = HARMONIC_ONSETS.get(word)
if parent:
return doubling_parameter(parent, digits)
return solve_multiplier(word, RR(1), digits)
def doubling_parameter(word, digits=DIGITS):
return solve_multiplier(word, RR(-1), digits)
def partition_matrix(word, digits=DIGITS):
p = PERIOD[word]
r = superstable_parameter(word, digits)
c = RR(1) / RR(2)
orbit = [c]
y = c
for _ in range(1, p):
y = f_value(r, y)
orbit.append(y)
points = sorted([RR(0), RR(1)] + orbit)
unique = []
tolerance = RR(10) ** (-(int(digits) + 20))
for point in points:
if not unique or abs(point - unique[-1]) > tolerance:
unique.append(point)
n = len(unique) - 1
matrix = [[0 for _ in range(n)] for _ in range(n)]
for i in range(n):
a, b = unique[i], unique[i + 1]
image_a = f_value(r, a)
image_b = f_value(r, b)
lo = min(image_a, image_b)
hi = max(image_a, image_b)
for j in range(n):
midpoint = (unique[j] + unique[j + 1]) / 2
if lo < midpoint < hi:
matrix[i][j] = 1
return matrix
def matmul(left, right):
n = len(left)
return [[sum(left[i][k] * right[k][j] for k in range(n))
for j in range(n)] for i in range(n)]
def charpoly_coeffs(matrix):
"""Characteristic polynomial coefficients, high degree first."""
n = len(matrix)
b = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
coeffs = [1]
for k in range(1, n + 1):
ab = matmul(matrix, b)
trace = sum(ab[i][i] for i in range(n))
if trace % k:
raise ArithmeticError("nonintegral characteristic coefficient")
c = -trace // k
coeffs.append(c)
b = [[ab[i][j] + (c if i == j else 0) for j in range(n)]
for i in range(n)]
return coeffs
def poly_eval(coeffs, value):
out = RR(0)
for coefficient in coeffs:
out = out * value + RR(coefficient)
return out
def perron_root(coeffs, digits=DIGITS):
lo = RR(1) + RR(10) ** (-(int(digits) + 10))
hi = RR(2)
flo = poly_eval(coeffs, lo)
fhi = poly_eval(coeffs, hi)
if flo * fhi > 0:
raise ArithmeticError("Perron root is not bracketed")
for _ in range(6 * int(digits)):
mid = (lo + hi) / 2
fmid = poly_eval(coeffs, mid)
if flo * fmid <= 0:
hi = mid
fhi = fmid
else:
lo = mid
flo = fmid
return (lo + hi) / 2
def entropy_parameter(word, digits=DIGITS):
key = ("entropy", word, WORKING_BITS)
if key in _CACHE:
return _CACHE[key]
matrix = partition_matrix(word, digits)
root = perron_root(charpoly_coeffs(matrix), digits)
value = root.log()
_CACHE[key] = value
return value
def r_value(word, expression, digits=DIGITS):
if expression == "onset":
return onset_parameter(word, digits)
if expression == "superstable":
return superstable_parameter(word, digits)
if expression == "doubling":
return doubling_parameter(word, digits)
if expression == "entropy":
return entropy_parameter(word, digits)
raise ValueError("unknown expression %r" % (expression,))
def endpoint_comment(word, expression, value):
parts = ["period $%d$" % PERIOD[word]]
special = SPECIAL_ENDPOINT_COMMENTS.get((word, expression))
if special:
parts.append(special)
parent = HARMONIC_ONSETS.get(word)
if expression == "onset" and parent:
parts.append(
r"$r_{\mathrm{on}}$ is the first period-doubling point of "
r"$\mathtt{%s}$" % parent
)
return "; ".join(parts) + "."
def normalisation_comment(word, base):
"""The remark on a `c` row: which `r` row it converts."""
symbol = {"onset": r"r_{\mathrm{on}}",
"superstable": "r_W",
"doubling": r"r_{\mathrm{pd}}"}[base]
return (r"period $%d$; the same window's $%s$ in the $z\mapsto z^2+c$ "
r"normalisation CITE{formula-conversion}."
% (PERIOD[word], symbol))
def entropy_comment(word):
polynomial = ENTROPY_POLYNOMIALS[word]
pieces = [
"period $%d$" % PERIOD[word],
r"if $\lambda$ is the Perron root then $\lambda$ satisfies $%s$" % polynomial,
r"$h_{\mathrm{top}}=\log\lambda$",
]
if word == "RLC":
pieces.append(
r"$\lambda$ is HREF{Golden_ratio#phi}[the golden ratio $\varphi$]")
if word == "RLRRRC":
pieces.append(r"$h_{\mathrm{top}}$ is half the entropy of $\mathtt{RLC}$")
if word in HARMONIC_ONSETS:
pieces.append(r"the entropy equals that of $\mathtt{%s}$" %
HARMONIC_ONSETS[word])
return "; ".join(pieces) + "."
class LogisticPeriodicWindows(numberdb.Generator):
table = os.environ.get("NUMBERDB_TABLE") or "T169"
parameters = ("W", "quantity")
type = "R"
digits = DIGITS
rigour = "heuristic (agreement-checked)"
def enumerate(self):
for word in WINDOWS:
for quantity in QUANTITIES:
yield {"W": word, "quantity": quantity}
def value(self, params, digits):
word = str(params["W"])
quantity = str(params["quantity"])
if word not in WINDOWS:
raise ValueError("W must be one of the listed kneading words")
if quantity not in QUANTITIES:
raise ValueError("quantity must be one of %s" %
", ".join(QUANTITIES))
base = BASE_OF.get(quantity, quantity)
value = r_value(word, base, digits)
if quantity in BASE_OF:
value = c_parameter(value)
comment = normalisation_comment(word, base)
elif quantity == "entropy":
comment = entropy_comment(word)
else:
comment = endpoint_comment(word, quantity, value)
exact = EXACT.get((word, quantity))
if exact is None:
return {"number": decimal(value, digits), "comment": comment}
if abs(RR(QQ(exact)) - RR(value)) >= RR(2) ** (-(WORKING_BITS // 2)):
raise ValueError(
"%s in %s is written exactly as %s, and the solve does not "
"agree" % (word, quantity, exact))
return {"number": exact, "comment": comment}
if __name__ == "__main__":
if os.environ.get("NUMBERDB_KEY_FROM_STDIN") == "1":
os.environ["NUMBERDB_API_KEY"] = sys.stdin.read().strip()
generator = LogisticPeriodicWindows()
if "--publish" in sys.argv or os.environ.get("NUMBERDB_PUBLISH") == "1":
print(generator.publish(
message="logistic-map periodic windows through period 8"))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)