rates.py
Checks: RATES: Walls sit at the preimages of the rationals criterion; RATES: The commitment bound theorem; RATES: Rate forcing theorem; RATES: The unit-index law pattern; RATES: The unit criterion criterion; RATES: The ladder law theorem.
Cited on: RATES.
To run it: save rates.py; standard library only.
"""
rates.py -- what a digit-by-digit reader pays at the continued-fraction
numeration, where the scale of a digit varies from stream to stream.
QUESTION. The continued-fraction numeration writes an irrational y as
its partial quotients [a_0; a_1, a_2, ...]. Its depth-n TILE, the
irrationals sharing the first n quotients, is the open interval between
the convergent p_(n-1)/q_(n-1) and (p_(n-1) + p_(n-2))/(q_(n-1) +
q_(n-2)), of length 1/(q_(n-1)(q_(n-1) + q_(n-2))). Unlike a positional
numeration, the SCALE of a tile, ln(1/length), is not a function of
its depth: one quotient buys 2 ln(phi) of scale on the all-ones stream
and far more where quotients are large. A reader of a monotone map f
sees the input's depth-n tile I_n and commits the longest quotient
string shared by every point of f(I_n); e(n) is its length, and the
LAG is n - e(n). What decides the lag, which maps have bounded lag, and
where does emission stop for good?
THE ARGUMENT (written before this script).
(1) WALLS ARE THE RATIONALS' PREIMAGES. The tiles are open intervals
whose endpoints are rationals, and every rational v = [b_0; ...,
b_k] (b_k >= 2 when k >= 1) separates two adjacent tiles at depth
k + 1, the streams near v agreeing on b_0..b_(k-1) and differing
at index k. Take f continuous and strictly monotone near an
irrational y. If f(y) is irrational, f(I_n) shrinks to a point
interior to each of its tiles, so emission grows without bound.
If f(y) = v is rational, every f(I_n) holds v in its interior,
and emission freezes at exactly k. So f walls exactly on f^-1(Q).
An integer Mobius map sends irrationals to irrationals: no walls.
y -> y^2 walls at the square roots of the non-square positive
rationals, and y -> sqrt(y) nowhere.
(2) THE COMMITMENT BOUND. Let z be irrational with a_j <= A for every
j >= 1, J an interval holding z, and C_e the deepest tile holding
J. The DEFICIT scale(J) - scale(C_e) is at most 3 ln(A + 2) +
ln 4, whatever J is. Proof: every rational p/q has |z - p/q| >=
1/((A + 2) q^2) (at a convergent from |z - p_n/q_n| =
1/(q_n(z_(n+1) q_n + q_(n-1))), elsewhere by Legendre's
criterion); a tile C_n has endpoints of denominator at most
2 q_(n-1) and length at most 1/q_(n-1)^2, so z sits at distance
at least |C_n|/(4(A + 2)) from its boundary. J fails to fit
C_(e+1), so |J| >= |C_(e+1)|/(4(A + 2)), and |C_e|/|C_(e+1)| <=
(A + 1)(A + 2).
(3) RATE FORCING. Let y have an eventually periodic tail of period
l and period matrix P, the product of [[a, 1], [1, 0]] over one
period, with dominant eigenvalue eta. Then ln q_n = n ln(eta)/l
+ O(1), so scale(I_n) = n r(y) + O(1) with the RATE r(y) =
2 ln(eta)/l. If f is C^1 and monotone with f'(y) != 0 and z =
f(y) is quadratic, scale(f(I_n)) = scale(I_n) + O(1), and (2)
gives e(n) r(z) = n r(y) + O(1): e(n) = rho n + O(1) with rho =
r(y)/r(z). The lag is bounded iff the rates match, and grows
linearly otherwise, with the sign of 1 - rho.
(4) THE CONDUCTOR LAW. eta is the fundamental unit of the multiplier
ring O of the lattice Z + Zy, which is eps^i, eps the fundamental
unit of the field and i = [O_max* : O*] the unit index, O_max the
ring of integers of Q(y); so r(y) = 2 i ln(eps)/l, and within one
field rho = (i_in l_out)/(i_out l_in) is rational, the logarithm
of the unit cancelling. That
the period generates the stabilizer is classical; this script
tests eta = eps^i with eps and i computed without the period.
(5) THE UNIT CRITERION. Every integer Mobius map with det = +-1 has
bounded lag at every irrational stream. It is a word in x + 1,
-x and 1/x, each a bounded-lookahead rewriting of the quotient
string: x + 1 raises a_0; 1/x moves the string by one place
(for x > 0, and 1/x = -(1/(-x)) for x < 0); and -[a_0; a_1, a_2,
...] = [-a_0 - 1; 1, a_1 - 1, a_2, ...] when a_1 >= 2, =
[-a_0 - 1; a_2 + 1, a_3, ...] when a_1 = 1. The lags of a fixed
word add. With (3) this makes r invariant under GL(2, Z). If
|det| = N >= 2, M = U diag(1, N) V (Smith form, U and V
unimodular) has lag slope 1 - r(phi)/r(phi/N) at y = V^-1 phi.
r(phi) = 2 ln(phi) is the least rate, attained only by an
all-ones tail (Perron-Frobenius: a period matrix dominates the
all-ones one entrywise). An all-ones tail would make phi/N a
unimodular image of phi, with the same multiplier ring Z[phi],
but the ring of Z + Z(phi/N) omits phi. So the lag grows: phi
is a universal witness.
(6) THE LADDER LAW. At a Stern-Brocot vertex v with parents l and
r, the straddle S_k(v) runs between the k-th mediants (p_l +
k p_v)/(q_l + k q_v) and (p_r + k p_v)/(q_r + k q_v), and holds v
in its interior at every k. For a point z on v's side of the far
parent (q_far its denominator), z lies in S_k(v) for exactly
K(v, z) = floor(1/(q_v^2 |z - v|) - q_far/q_v) values k >= 1.
At a convergent p_n/q_n of z the far parent is p_(n-1)/q_(n-1),
and K = a_(n+1). At any other vertex of z's path Legendre gives
|z - v| >= 1/(2 q_v^2), so K <= 1. At z = v the chain is
infinite. Adding the straddles to the tiles (the cover of the
redundant digit M) extends exactly the terminated ladders: it
cures walls and adds no rank elsewhere, so it cannot cure a
rate mismatch.
HAND-ATTACK (on paper, before any engine code).
- e counts a_0. v = 2 = [2] separates the depth-1 tiles (1, 2) and
(2, 3), so the frozen emission at a wall on 2 is 0; 3/2 = [1; 2]
freezes at 1, 7/3 = [2; 3] at 1, 17/12 = [1; 2, 2, 2] at 3.
- The deficit constant, re-derived: dist(z, boundary C_n) >=
1/((A + 2)(2 q)^2) >= |C_n|/(4(A + 2)), so the bound is
ln((A + 1)(A + 2)) + ln(4(A + 2)) <= 3 ln(A + 2) + ln 4.
- rho by hand, from (i, l): phi (1, 1); 2 phi = [3; (4)] has eta =
2 + sqrt5 = phi^3, (3, 1); 3 phi = [4; (1, 5)], eta = phi^4, (4,
2); phi/2 = [0; 1, (4)], (3, 1); 7 phi = [11; (3, 15)], eta =
phi^8, (8, 2); sqrt2 (1, 1), 2 sqrt2 (1, 1) [Pell], 5 sqrt2 =
[7; (14)], (3, 1); sqrt3 = [1; (1, 2)], (1, 2), 3 sqrt3 (3, 2);
1 + sqrt3 (1, 2), (1 + sqrt3)^2 = 4 + 2 sqrt3 = [7; (2, 6)], (2,
2); 7 + 4 sqrt3 = [13; (1, 12)], lattice ring Z[4 sqrt3], (2, 2),
and its root 2 + sqrt3 (1, 2). So rho = 1/3 at (2x, phi), 1/2 at
(3x, phi), 1/3 at (x/2, phi), 1 at (2x, sqrt2), 1/3 at (5x,
sqrt2), 1/3 at (3x, sqrt3), 1 at (x^2, phi), 1/2 at (x^2, 1 +
sqrt3), 2 at (sqrt x, 7 + 4 sqrt3), 3 at (x/2, 2 phi), 1/4 at
(7x, phi).
- A wall's straddle index: at v = 2 the chain's half-widths are
1/(k + 1) and 1/k, so k(n) ~ 1/|f(I_n)| and ln k(n) grows at the
input's rate: 2 ln(1 + sqrt2) = 1.7627 per quotient at sqrt2.
- Counting in Stern-Brocot steps instead of quotients (a quotient a
is a steps), (2x, phi) emits 4 steps per 3 input steps, since 2 phi
has quotients 4: the lag's SIGN depends on the unit of count, the
mismatch does not.
PREDICTIONS (frozen before any engine code).
P-0 Positive controls, read first: 1/x at phi has lag <= 1 at every
depth; 2x at phi has e(150)/150 within 0.02 of 1/3; y^2 at phi
emits at least 100 quotients from 120.
P-W y^2 freezes at input depths 20, 60 and 120 at sqrt2, sqrt3,
sqrt6 (at 0), sqrt(3/2), sqrt(7/3) (at 1) and sqrt(17/12) (at
3), each the index of its image's last quotient; y^2 emits
without bound at phi and 1 + sqrt2, sqrt(y) at 3 + sqrt2.
P-D At (2x, phi), (2x, sqrt2), (y^2, phi) and (7x, phi), with A =
4, 4, 1, 15, the deficit lies in [0, 3 ln(A + 2) + ln 4] at
every depth to 240 from the first emission, and so does every
interval of a random battery around sqrt7 and 7 phi. At the
wall (y^2, sqrt2) the deficit grows with slope 1.7627 +- 5%.
P-F At the eleven rows above, rho from the lattice data equals the
hand value, the residual |e(n) - rho n| over depths 121..240
stays within its sup over 1..120 plus 2, and at the mismatched
rows the final lag has the sign of 1 - rho.
P-C eta = eps^i exactly at every scanned quadratic y (at least 50,
over fields of discriminant up to 100), with eps found by
searching x^2 - d y^2 = +-1 (or +-4) and i the least power of
eps in the multiplier ring.
P-U x + 1, x - 1, 1/x, -x, -1/x, (2x + 1)/(x + 1) and the
sign-crossing (x - 4)/(x - 3) have lag bounded by 4 and
non-growing at phi, sqrt2, sqrt3 and a fixed pseudo-random
string of quotients 1..5, and keep (eta, l) at every quadratic
input. N phi has no all-ones tail and a rate above 2 ln(phi) at
every N = 2..40.
P-K The closed form for K equals enumeration at every vertex of the
scanned paths; K = a_(n+1) at every convergent and K <= 1 at
every other path vertex. At the walls of y^2 at sqrt2, sqrt3 and
sqrt(3/2) the straddle reader never freezes, its ln k slope
within 3% of the input's rate. In Stern-Brocot steps (2x, phi)
emits at 4/3 and (x/2, 2 phi) at 3/4, each within 4%.
FINDINGS (entered after the run; prints copied from it).
- The first run passed 16 of 18. The two ladder checks failed on a
bound in the script, not in the law: the path walk ran past the last
convergent computed, so later convergents were
counted as "other" vertices and the last one had no a_(n+1). The
walk now stops at the last convergent whose next quotient is known.
The closed form had agreed with enumeration throughout.
- The hand values listed (1, 1) for 2 sqrt2. The script prints (2, 2):
2 sqrt2 = [2; (1, 4)] with eta = 3 + 2 sqrt2 = (1 + sqrt2)^2, so the
Pell row's rho = 1 comes from the index and the period doubling
together. rho was right; the pair was not.
- P-0 holds: 1/x at phi has lag -1 at depths 1 to 150, 2x at phi emits 50
quotients from 150, and y^2 at phi emits 118 from 120.
- P-W holds: y^2 freezes at [0, 0, 0] at sqrt2, sqrt3 and sqrt6, at
[1, 1, 1] at sqrt(3/2) and sqrt(7/3), and at [3, 3, 3] at
sqrt(17/12). It emits [18, 58, 118] at phi and at 1 + sqrt2, and
sqrt(y) emits [10, 38, 86] at 3 + sqrt2.
- P-D holds: deficits lie in [0, 2.015] at (2x, phi), [0.405, 2.199]
at (2x, sqrt2), [0.639, 0.831] at (y^2, phi) and [0.754, 5.032] at
(7x, phi), against bounds 6.762, 6.762, 4.682 and 9.886. All 296
random intervals that commit a tile (of 300) obey the bound, the
largest excess -3.235. At the wall the deficit's slope is
1.7627, and 2 ln(1 + sqrt2) = 1.7627.
- P-F holds at all eleven rows. The residual sup is the same over
both halves at every row (at most 2.00, at (y^2, phi)). Final lags
at 240: 160, 120, 159, 0, 160, 160, 2, 121, -239, -479 and 180.
- P-C holds: eta = eps^i at 130 distinct quadratic points over 15 fields,
conductors up to 20.
- P-U holds: the largest lag of any unimodular map at any input is 1
(-x at phi and sqrt3); each keeps (eta, l). N phi's rate is 2.000
to 9.000 times r(phi) over N = 2..40.
- P-K holds: the closed form equals enumeration at 531 convergents
and 1138 other path vertices, K = a_(n+1) at every convergent and
K <= 1 elsewhere. The straddle reader's ln k slope is 1.7627,
1.3170 and 2.2924 at sqrt2, sqrt3 and sqrt(3/2), each equal to the
input's rate to four places. In Stern-Brocot steps (2x, phi) emits
at 1.3292 and (x/2, 2 phi) at 0.7497.
SETTLED ON AUDIT (the argument above stays as frozen; these replace it).
- (5) holds for matrices with coprime entries. A common factor g
multiplies det by g^2 and leaves the map alone (2y/2 is the
identity, det 4, lag 0), and the Smith form diag(1, N) needs the
entries coprime.
- (6) counts straddles for z != v inside v's tile (l, r); for z
outside it the floor goes negative while the count is 0.
- (6) "the cover of the redundant digit M": what is checked is that
S_1 of the root is (1/2, 2), the image of [0, inf] under
y -> (2y + 1)/(y + 2).
- (5) bounds the lag above by the word's lookahead; below, M^-1 is a
word too, and a depth-n tile mapped into a depth-e tile forces
e <= n + c. The check now bounds min(lag) as well: every unimodular
row's lags lie in [-2, 1].
- (6) holds for irrational z > 0: at a rational z on a mediant the
floor over-counts by one (v = 1, z = 2/3 gives 2, but z lies only
in S_1). The parents of v are the two Farey neighbours v is the
mediant of; at an integer vertex one has an equal denominator.
- (4) is a unit-index law, not a conductor law: i does not fix the
conductor (eps = 8 + 3 sqrt7 lies in Z[3 sqrt7], so i = 1 at
conductor 3).
- (x - 4)/(x - 3) crosses no sign at these inputs, all in (1, 2): it
exercises a negative denominator only.
- (6)'s "cannot cure a rate mismatch": no reader commits past f(I_n),
of scale n r(y) + O(1); the cover is run only at the walls.
- [F]'s rho is r(y)/r(z) from the period eigenvalues (rate), not
from (i, l), so rate forcing's check does not lean on (4).
RUN RECORD: 18/18 checks, 0.18 s, peak commit 8 MB under a memory guard.
Later the wall indices are read off y^2's own quotients; the
unimodular lags are asserted in [-2, 1], the residual's two sups equal
and at most 2, the lag's sign strictly where rho != 1, P-C's 50 points
and S_1 of the root checked: 20/20, 0.2 s, peak commit 8 MB. Last,
[F]'s rho from the period eigenvalues, the check names scoped to their
loops and P-C's fixed count dropped: 20/20, 0.2 s, peak commit 8 MB
(rho times 1.001 fails the rho check, in a scratch copy).
"""
import math
import random
CHECKS = []
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
tail = f" ({detail})" if detail else ""
print(f" {'PASS' if ok else 'FAIL'} {name}{tail}")
# ---- quadratic numbers (u + v sqrt d)/w, d not a square, w > 0 ----
def qn(u, v, w, d):
if w < 0:
u, v, w = -u, -v, -w
g = math.gcd(math.gcd(u, v), w)
return (u // g, v // g, w // g, d)
def qfloor(x):
u, v, w, d = x
s = math.isqrt(v * v * d)
top = u + s if v > 0 else (u - s - 1 if v < 0 else u)
return top // w
def qsub(x, a):
u, v, w, d = x
return qn(u - a * w, v, w, d)
def qrecip(x):
u, v, w, d = x
return qn(w * u, -w * v, u * u - v * v * d, d)
def qmob(x, m):
a, b, c, e = m
u, v, w, d = x
nu, nv, du, dv = a * u + b * w, a * v, c * u + e * w, c * v
return qn(nu * du - nv * dv * d, nv * du - nu * dv,
du * du - dv * dv * d, d)
def qsq(x):
u, v, w, d = x
return qn(u * u + v * v * d, 2 * u * v, w * w, d)
def qcf(x, n):
out = []
for _ in range(n):
a = qfloor(x)
out.append(a)
x = qrecip(qsub(x, a))
return out
def qperiod(x):
seen, quots = {}, []
while x not in seen:
seen[x] = len(quots)
a = qfloor(x)
quots.append(a)
x = qrecip(qsub(x, a))
s = seen[x]
return quots[:s], quots[s:]
def qsign(x, p, q):
"""Sign of x - p/q, exactly (q > 0)."""
u, v, w, d = x
a, b = q * u - p * w, q * v
if a >= 0 and b >= 0:
return 1 if a or b else 0
if a <= 0 and b <= 0:
return -1
t = a * a - b * b * d
return (1 if t > 0 else -1) if a > 0 else (1 if t < 0 else -1)
# ---- lattice data: eta from the period, eps and i without it ----
def squarefree(n):
s, k = 1, 2
while k * k <= n:
while n % (k * k) == 0:
n //= k * k
s *= k
k += 1
return n, s
def period_eta(per):
"""(tr, det) of the period matrix."""
a0, b0, c0, d0 = 1, 0, 0, 1
for a in per:
a0, b0, c0, d0 = a0 * a + b0, a0, c0 * a + d0, c0
return a0 + d0, a0 * d0 - b0 * c0
def eta_exact(per, d0):
"""eta = (tr + s sqrt d0)/2 as the integer pair (tr, s)."""
tr, det = period_eta(per)
disc = tr * tr - 4 * det
assert disc % d0 == 0
s = math.isqrt(disc // d0)
assert s * s * d0 == disc
return tr, s
def fund_unit(d0):
"""Least unit > 1 of the maximal order, as (X, Y): (X + Y sqrt d0)/2."""
y = 1
while True:
for sgn in (-1, 1):
if d0 % 4 == 1:
t = d0 * y * y + 4 * sgn
x = math.isqrt(t) if t > 0 else -1
if x >= 0 and x * x == t:
return x, y
else:
t = d0 * y * y + sgn
x = math.isqrt(t) if t > 0 else -1
if x >= 0 and x * x == t:
return 2 * x, 2 * y
y += 1
def hmul(a, b, d0):
"""(X1 + Y1 s)/2 * (X2 + Y2 s)/2 with s = sqrt d0, halves kept."""
x = a[0] * b[0] + a[1] * b[1] * d0
y = a[0] * b[1] + a[1] * b[0]
assert x % 2 == 0 and y % 2 == 0
return x // 2, y // 2
def in_order(el, d0, f):
"""Is (X + Y sqrt d0)/2 in Z + f O_max?"""
x, y = el
if d0 % 4 == 1:
return (x - y) % 2 == 0 and y % f == 0
return x % 2 == 0 and y % 2 == 0 and (y // 2) % f == 0
def lattice(x):
"""(d0, per, f, eta, eps^i, i) for a quadratic x = (u + v sqrt d)/w."""
u, v, w, d = x
d0, s = squarefree(d)
aa, bb, cc = w * w, -2 * u * w, u * u - v * v * d
g = math.gcd(math.gcd(aa, bb), cc)
disc = (bb * bb - 4 * aa * cc) // (g * g)
dk = d0 if d0 % 4 == 1 else 4 * d0
f = math.isqrt(disc // dk)
assert f * f * dk == disc
_, per = qperiod(x)
eta = eta_exact(per, d0)
eps = fund_unit(d0)
pw, i = eps, 1
while not in_order(pw, d0, f):
pw, i = hmul(pw, eps, d0), i + 1
return d0, per, f, eta, pw, i
def rate(per, d0):
tr, s = eta_exact(per, d0)
return 2 * math.log((tr + s * math.sqrt(d0)) / 2) / len(per)
# ---- the interval reader ----
def conv(qs):
p, q, pp, qq = 1, 0, 0, 1
for a in qs:
p, q, pp, qq = a * p + pp, a * q + qq, p, q
return p, q, pp, qq
def tile(qs):
p, q, pp, qq = conv(qs)
e1, e2 = (p, q), (p + pp, q + qq)
return (e1, e2) if e1[0] * e2[1] < e2[0] * e1[1] else (e2, e1)
def frac_order(a, b):
return (a, b) if a[0] * b[1] < b[0] * a[1] else (b, a)
def mob(m):
a, b, c, e = m
def f(lo, hi):
out = []
for n, dd in (lo, hi):
nn, nd = a * n + b * dd, c * n + e * dd
out.append((-nn, -nd) if nd < 0 else (nn, nd))
return frac_order(*out)
return f
def square(lo, hi):
return frac_order((lo[0] ** 2, lo[1] ** 2), (hi[0] ** 2, hi[1] ** 2))
SCALE = 10 ** 400
def root(lo, hi):
a = math.isqrt(lo[0] * SCALE * SCALE // lo[1])
b = math.isqrt(-(-hi[0] * SCALE * SCALE // hi[1])) + 1
return (a, SCALE), (b, SCALE)
def emit(lo, hi):
(ln, ld), (hn, hd) = lo, hi
out = []
while True:
a = ln // ld
if hn > (a + 1) * hd:
return out
out.append(a)
if ln == a * ld:
return out
ln, ld, hn, hd = hd, hn - a * hd, ld, ln - a * ld
def scale_of_tile(qs):
_, q, _, qq = conv(qs)
return math.log(q * (q + qq))
def scale_of(lo, hi):
return math.log(lo[1] * hi[1]) - math.log(hi[0] * lo[1] - lo[0] * hi[1])
def reader(stream, f, n):
lo, hi = f(*tile(stream[:n]))
return lo, hi, emit(lo, hi)
PHI = (1, 1, 2, 5)
SQ2, SQ3 = (0, 1, 1, 2), (0, 1, 1, 3)
def lcg_stream(n, seed=12345):
out, s = [1], seed
for _ in range(n - 1):
s = (1103515245 * s + 12345) % 2 ** 31
out.append(1 + s % 5)
return out
def section_0():
print("\n[0] positive controls")
st = qcf(PHI, 160)
lags = [n - len(reader(st, mob((0, 1, 1, 0)), n)[2])
for n in range(1, 151)]
e150 = len(reader(st, mob((2, 0, 0, 1)), 150)[2])
e120 = len(reader(st, square, 120)[2])
print(f" 1/x at phi: lag range [{min(lags)}, {max(lags)}]; "
f"2x at phi: e(150) = {e150}; y^2 at phi: e(120) = {e120}")
check("1/x at phi has lag <= 1 at depths 1..150", max(lags) <= 1)
check("2x at phi emits at slope 1/3", abs(e150 / 150 - 1 / 3) < 0.02)
check("y^2 at phi emits (the wall detector sees no wall)", e120 >= 100)
def section_w():
print("\n[W] walls sit at the preimages of the rationals")
rows = [("sqrt2", SQ2), ("sqrt3", SQ3), ("sqrt6", (0, 1, 1, 6)),
("sqrt(3/2)", (0, 1, 2, 6)), ("sqrt(7/3)", (0, 1, 3, 21)),
("sqrt(17/12)", (0, 1, 6, 51))]
ok = True
for name, y in rows:
num, den, k = y[3], y[2] ** 2, -1 # y^2 = d/w^2, its quotients
while den:
num, den, k = den, num % den, k + 1
st = qcf(y, 130)
es = [len(reader(st, square, n)[2]) for n in (20, 60, 120)]
print(f" y^2 at {name:12s}: emission {es}, predicted {k}")
ok &= es == [k] * 3
check("y^2 freezes at each wall at its image's last-quotient index", ok)
grow = []
for name, y, f in [("y^2 at phi", PHI, square),
("y^2 at 1 + sqrt2", (1, 1, 1, 2), square),
("sqrt(y) at 3 + sqrt2", (3, 1, 1, 2), root)]:
st = qcf(y, 130)
es = [len(reader(st, f, n)[2]) for n in (20, 60, 120)]
print(f" {name:22s}: emission {es}")
grow.append(es[0] < es[1] < es[2])
check("off the rationals emission grows over depths 20, 60, 120",
all(grow))
def deficits(stream, f, n_max, whole=False):
"""Before the first emission no tile is committed; with whole=True
the committed set is the whole line, of scale 0."""
out = []
for n in range(1, n_max + 1):
lo, hi, e = reader(stream, f, n)
if e:
out.append(scale_of(lo, hi) - scale_of_tile(e))
else:
out.append(scale_of(lo, hi) if whole else None)
return out
def section_d():
print("\n[D] the commitment bound")
rows = [("2x at phi", PHI, mob((2, 0, 0, 1)), 4),
("2x at sqrt2", SQ2, mob((2, 0, 0, 1)), 4),
("y^2 at phi", PHI, square, 1),
("7x at phi", PHI, mob((7, 0, 0, 1)), 15)]
ok = True
for name, y, f, a in rows:
ds = [x for x in deficits(qcf(y, 250), f, 240) if x is not None]
bound = 3 * math.log(a + 2) + math.log(4)
print(f" {name:12s}: deficit in [{min(ds):.3f}, {max(ds):.3f}],"
f" bound {bound:.3f} (A = {a})")
ok &= min(ds) >= -1e-9 and max(ds) <= bound
check("the deficit stays in [0, 3 ln(A + 2) + ln 4] at every row", ok)
rng = random.Random(7)
worst, count = -math.inf, 0
for z, a in [((0, 1, 1, 7), 4), ((7, 7, 2, 5), 15)]:
qs = qcf(z, 40)
for _ in range(150):
j = rng.randrange(2, 30)
p1, q1, _, _ = conv(qs[:j])
p2, q2, _, _ = conv(qs[:j + 1])
lo, hi = frac_order((p1, q1), (p2, q2))
den = q1 * q1 * rng.randrange(1, 50)
lo = (lo[0] * den - lo[1] * rng.randrange(0, 5), lo[1] * den)
hi = (hi[0] * den + hi[1] * rng.randrange(0, 5), hi[1] * den)
e = emit(lo, hi)
if e:
gap = scale_of(lo, hi) - scale_of_tile(e)
worst = max(worst, gap - 3 * math.log(a + 2) - math.log(4))
count += 1
print(f" random battery: {count} of 300 intervals around sqrt7 and"
f" 7 phi commit a tile (the rest nothing), largest excess over"
f" the bound {worst:.3f}")
check("every random interval obeys the bound", worst <= 0 and count > 250)
ds = deficits(qcf(SQ2, 250), square, 240, whole=True)
slope = (ds[239] - ds[119]) / 120
print(f" wall y^2 at sqrt2: deficit slope {slope:.4f} against "
f"2 ln(1 + sqrt2) = {2 * math.log(1 + math.sqrt(2)):.4f}")
check("at the wall the deficit grows at the input's rate",
abs(slope / (2 * math.log(1 + math.sqrt(2))) - 1) < 0.05)
def section_f():
print("\n[F] rate forcing")
rows = [("2x", PHI, (2, 0, 0, 1), 1 / 3), ("3x", PHI, (3, 0, 0, 1), 1 / 2),
("x/2", PHI, (1, 0, 0, 2), 1 / 3), ("2x", SQ2, (2, 0, 0, 1), 1),
("5x", SQ2, (5, 0, 0, 1), 1 / 3), ("3x", SQ3, (3, 0, 0, 1), 1 / 3),
("y^2", PHI, "sq", 1), ("y^2", (1, 1, 1, 3), "sq", 1 / 2),
("sqrt y", (7, 4, 1, 3), "rt", 2),
("x/2", (1, 1, 1, 5), (1, 0, 0, 2), 3),
("7x", PHI, (7, 0, 0, 1), 1 / 4)]
all_ok, same = True, True
for name, y, m, hand in rows:
if m == "sq":
f, z = square, qsq(y)
elif m == "rt":
f, z = root, (2, 1, 1, 3)
else:
f, z = mob(m), qmob(y, m)
ly, lz = lattice(y), lattice(z)
rho = rate(ly[1], ly[0]) / rate(lz[1], lz[0])
st = qcf(y, 250)
es = [len(reader(st, f, n)[2]) for n in range(1, 241)]
res = [abs(es[n - 1] - rho * n) for n in range(1, 241)]
first, second = max(res[:120]), max(res[120:])
lag = 240 - es[-1]
sign_ok = (abs(rho - 1) < 1e-9
or (lag > 0 if rho < 1 else lag < 0))
ok = abs(rho - hand) < 1e-9 and second <= first + 2 and sign_ok
all_ok &= ok
same &= abs(second - first) < 1e-9 and first <= 2
print(f" ({name}, y = {y[:3]}/{y[3]}): (i, l) {ly[5], len(ly[1])} "
f"-> {lz[5], len(lz[1])}, rho {rho:.4f}, residual sup "
f"{first:.2f} then {second:.2f}, lag(240) {lag}")
check("rho matches the hand value, e(n) - rho n stays flat, and the "
"lag's sign is 1 - rho, at all eleven rows", all_ok)
check("the residual's sup is the same over both halves, at most 2",
same)
def section_c():
print("\n[C] the unit-index law: eta = eps^i")
ys = []
for d in (2, 3, 5, 6, 7, 10, 11, 13, 14, 15, 17, 19, 21, 22, 23):
for (u, v, w) in [(0, 1, 1), (1, 1, 1), (1, 1, 2), (0, 2, 1),
(1, 3, 1), (2, 1, 3), (0, 3, 2), (1, 2, 5)]:
ys.append(qn(u, v, w, d))
ys += [qn(n, n, 2, 5) for n in range(2, 13)]
ys = list(dict.fromkeys(ys)) # qn(2, 2, 2, 5) is d = 5's (1, 1, 1)
ok, fields, conductors = True, set(), set()
for y in ys:
d0, per, f, eta, pw, i = lattice(y)
ok &= eta == pw
fields.add(d0)
conductors.add(f)
print(f" {len(ys)} quadratic points over {len(fields)} fields, "
f"conductors up to {max(conductors)}")
check("the period's eigenvalue equals eps^i at every point", ok)
def section_u():
print("\n[U] the unit criterion")
maps = [("x + 1", (1, 1, 0, 1)), ("x - 1", (1, -1, 0, 1)),
("1/x", (0, 1, 1, 0)), ("-x", (-1, 0, 0, 1)),
("-1/x", (0, -1, 1, 0)), ("(2x + 1)/(x + 1)", (2, 1, 1, 1)),
("(x - 4)/(x - 3)", (1, -4, 1, -3))]
inputs = [("phi", qcf(PHI, 130), PHI), ("sqrt2", qcf(SQ2, 130), SQ2),
("sqrt3", qcf(SQ3, 130), SQ3), ("lcg", lcg_stream(130), None)]
bounded, kept = True, True
for name, m in maps:
worst = []
for iname, st, y in inputs:
lags = [n - len(reader(st, mob(m), n)[2]) for n in range(2, 121)]
h = len(lags) // 2
bounded &= (max(lags) <= 1 and min(lags) >= -2
and max(lags[h:]) <= max(lags[:h]) + 1)
worst.append((min(lags), max(lags)))
if y is not None:
a, b = lattice(y), lattice(qmob(y, m))
kept &= a[3] == b[3] and len(a[1]) == len(b[1])
print(f" {name:18s}: lag range at phi, sqrt2, sqrt3, lcg "
f"{worst}")
check("every unimodular map has bounded, non-growing lag", bounded)
check("every unimodular image keeps (eta, l)", kept)
ok, low = True, []
for n in range(2, 41):
_, per = qperiod(qn(n, n, 2, 5))
r = rate(per, 5)
ok &= set(per) != {1} and r > 2 * math.log((1 + 5 ** 0.5) / 2)
low.append(r / (2 * math.log((1 + 5 ** 0.5) / 2)))
print(f" N phi, N = 2..40: rate over r(phi) from {min(low):.3f}"
f" to {max(low):.3f}")
check("N phi has no all-ones tail and a rate above r(phi)", ok)
def sb_parents(p, q):
l, r = (0, 1), (1, 0)
while True:
m = (l[0] + r[0], l[1] + r[1])
if m == (p, q):
return l, r
if p * m[1] < m[0] * q:
r = m
else:
l = m
def straddle_k(lo, hi, v, l, r):
(ln, ld), (hn, hd) = lo, hi
(pv, qv), (pl, ql), (pr, qr) = v, l, r
num, den = ln * ql - pl * ld, pv * ld - ln * qv
kl = -(-num // den) - 1 if num > 0 else 0
num, den = hd * pr - hn * qr, hn * qv - hd * pv
kr = -(-num // den) - 1 if num > 0 else 0
return min(kl, kr)
def section_k():
print("\n[K] the ladder law and the redundant cover")
l, r = sb_parents(1, 1)
s1 = sorted([(l[0] + 1, l[1] + 1), (r[0] + 1, r[1] + 1)],
key=lambda t: t[0] / t[1])
a, b, c, d = 2, 1, 1, 2 # y -> (2y + 1)/(y + 2)
image = sorted([(b, d), (a, c)], key=lambda t: t[0] / t[1])
print(f" S_1 of the root: {s1}; the image of 0 and inf: {image}")
check("S_1 of the root is the image of (0, inf) under (2y + 1)/(y + 2)",
s1 == image)
zs = [PHI, (1, 1, 1, 5), SQ2, SQ3, (0, 1, 1, 7), (3, 3, 2, 5),
(1, 1, 2, 13), (7, 7, 2, 5), (5, 1, 1, 6)]
agree, at_conv, elsewhere, n_conv, n_other = True, True, True, 0, 0
for z in zs:
qs = qcf(z, 60)
convs = {}
p, q, pp, qq = 1, 0, 0, 1
for j, a in enumerate(qs):
p, q, pp, qq = a * p + pp, a * q + qq, p, q
convs[(p, q)] = j
qmax = conv(qs[:-1])[1]
l, r = (0, 1), (1, 0)
while True:
v = (l[0] + r[0], l[1] + r[1])
if v[1] > qmax:
break
side = qsign(z, *v)
far = l if side < 0 else r
k = 0
while True:
m = (far[0] + (k + 1) * v[0], far[1] + (k + 1) * v[1])
if qsign(z, *m) != -side:
break
k += 1
u, vv, w, d = z
dist = qn(side * (u * v[1] - v[0] * w), side * vv * v[1],
w * v[1], d)
x = qrecip((dist[0] * v[1] ** 2, dist[1] * v[1] ** 2,
dist[2], d))
x = qn(x[0] * v[1] - far[1] * x[2], x[1] * v[1], x[2] * v[1], d)
agree &= qfloor(x) == k
if v in convs:
n_conv += 1
j = convs[v]
at_conv &= j + 1 < len(qs) and k == qs[j + 1]
else:
n_other += 1
elsewhere &= k <= 1
if side < 0:
r = v
else:
l = v
print(f" {len(zs)} points: {n_conv} convergents, {n_other} other "
f"path vertices")
check("the closed form for K equals enumeration at every vertex", agree)
check("K = a_(n+1) at every convergent", at_conv)
check("K <= 1 at every other path vertex", elsewhere)
ok = True
for name, y, v in [("sqrt2", SQ2, (2, 1)), ("sqrt3", SQ3, (3, 1)),
("sqrt(3/2)", (0, 1, 2, 6), (3, 2))]:
st = qcf(y, 130)
l, r = sb_parents(*v)
ks = [straddle_k(*square(*tile(st[:n])), v, l, r)
for n in (60, 120)]
_, per = qperiod(y)
rin = rate(per, y[3])
slope = (math.log(ks[1]) - math.log(ks[0])) / 60
print(f" straddle reader at y^2, {name:9s}: k(60), k(120) "
f"= {ks[0]:.3e}, {ks[1]:.3e}; ln k slope {slope:.4f}, "
f"input rate {rin:.4f}")
ok &= ks[0] > 1 and abs(slope / rin - 1) < 0.03
check("the straddle reader grows from depth 60 to 120 at a wall, at "
"the input's rate", ok)
rel = []
for y, m, want in [(PHI, (2, 0, 0, 1), 4 / 3),
((1, 1, 1, 5), (1, 0, 0, 2), 3 / 4)]:
st = qcf(y, 250)
e = reader(st, mob(m), 240)[2]
steps = sum(e) / sum(st[:240])
rel.append((steps, want))
print(f" Stern-Brocot steps emitted per step read: (2x, phi) "
f"{rel[0][0]:.4f}, (x/2, 2 phi) {rel[1][0]:.4f}")
check("in steps the lags change sign: 4/3 and 3/4",
all(abs(s / w - 1) < 0.04 for s, w in rel))
def main():
section_0()
section_w()
section_d()
section_f()
section_c()
section_u()
section_k()
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
raise SystemExit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()