poles.py
Checks: REDUNDANT: The pole units theorem; REDUNDANT: The selection estimate property; rule; REDUNDANT: The relative-gradient law property; rule.
Cited on: REDUNDANT.
To run it: save poles.py, with manystream.py, redundant.py beside it; standard library only.
"""
poles.py -- the pole units read by the many-stream theorem, and the
law that sets their lookahead: the square root, the reciprocal and the
divider, whose lookahead can be negative.
QUESTION. manystream.py reads every C^2 map of signed-digit streams
at the least L with b^L rho >= Lam |I|, Lam the sup over the root box of
the gradient's L1 norm, L = c + o the lookahead plus the LEAD. The
literature's on-line units with a pole -- the square root sqrt(s + x), the
reciprocal 1/(s + x) and the divider x/(s + y) -- are C^2 on the
root box once the pole sits outside it. Three questions. Does the
theorem read them when the lookahead c = L - o is NEGATIVE, the
reader emitting digits before it reads one? What sets the sign of c:
the pole, the dependency structure (a stream multiplied by a function
of the others), or one number of the map? And how do the literature's
published on-line delays sit against the floors the law gives?
THE READER is manystream.py's: streams of the digits D = {-a-, ...,
a+} in radix b, SLACK rho = a- + a+ + 1 - b >= 1, M+- = a+-/(b - 1),
w = M- + M+, |I| = a- + a+; output digits weigh b^(o - t); at
lookahead c the output digit t is committed after n = max(0, t + c)
digits of every stream. The POLE DISTANCE h > 0 places the pole:
s = M- + h, so s + x runs over [h, h + w] on the root box.
THE ARGUMENT (written before this script).
(1) EVERY L IN Z. For L < o the first o - L output digits are
committed on the root box alone. The level-t zone has length
(w - 1) b^(o - t) in absolute units, at least (w - 1) b^L at
every t <= o - L, and the root box's image is at most Lam w wide,
so under b^L rho >= Lam |I| no such digit is dead; below the law
the kill of manystream.py's necessity argument sits at a depth
n past which nothing changed. And a kill at a coarser level
implies one at the next: in units of the coarser level a zone
of length w - 1 holds a whole finer zone of length (w - 1)/b,
spaced 1/b, whenever
(w - 1)(b - 1) >= 1, which is rho >= 1. So the theorem holds
as stated at every L in Z, the lookahead of either sign.
(2) THE CLOSED FORMS. With Mh = max(M-, M+):
sqrt(s + x): Lam = 1/(2 sqrt h), range [sqrt h, sqrt(h + w)];
1/(s + x): Lam = 1/h^2, range [1/(h + w), 1/h];
x/(s + y): Lam = (h + Mh)/h^2, range [-M-/h, M+/h],
each partial's modulus largest at the pole corner y = -M-, the
divider's with |x| = Mh.
(3) THE LEAD IS THE RANGE'S SCALE. The root tile [-M- b^o, M+ b^o]
holds [lo, hi] iff b^o >= SIGMA = max(hi/M+, -lo/M-) (a+- >= 1),
so o = ceil(log_b SIGMA), while L* = ceil(log_b(Lam |I| / rho)).
Two ceilings differ by less than one more than their arguments,
so with the RELATIVE GRADIENT Gamma = Lam / SIGMA,
c* - log_b(Gamma |I| / rho) lies in (-1, 1):
THE RELATIVE-GRADIENT LAW. The map enters the lookahead through
Gamma alone, to within one digit.
(4) A MULTIPLIED STREAM. For f = x g(y), g > 0 on the root box with
sup G, the x-partial is g, so Lam >= G, and the range is
[-M- G, M+ G], so SIGMA = G and Gamma >= 1. With a+- <= b - 1,
rho <= b - 1 and |I|/rho = 1 + (b - 1)/rho >= 2, so
c* > log_b 2 - 1 > -1: c* >= 0. At slack 1, |I|/rho = b and
c* >= 1. The divider has Gamma = 1 + Mh/h, so c* < 2 +
log_b(1 + Mh/h), and c* <= 2 once Mh/h <= b - 1, which
h >= 1/(b - 1) grants since Mh <= 1.
(5) THE POLE UNITS RUN AHEAD. The root has Gamma =
M+/(2 sqrt(h (h + w))) and the reciprocal Gamma = M+/h, both
falling to 0 as h grows, so c* falls without bound: the reader
emits ever more digits before it reads one.
(6) THE STRUCTURE IS NOT THE VARIABLE. exp(-h y), one stream and no
product, has Lam = h e^(h M-) and SIGMA = e^(h M-)/M+, so
Gamma = h M+ RISES with h and c* climbs like log_b h. A stream
multiplied by a function of the others is one way to hold Gamma
up, not the only one.
THE DESIGN, frozen before the engine.
THE PAIRS: the 20 cells b = 2..5, 1 <= a+- <= b - 1, rho >= 1
(manystream.census), at h = j/(b - 1), j in {1, 2, 4, 8, 32, 128}.
THE CLOSED FORMS against manystream.Cell's integer loops (law_L,
lead) at every pair and unit; the root's decided by squaring.
THE SCANS at L* and L* - 1, L* the LAW'S least L (never raised to
the lead, unlike manystream.law_of): the reciprocal and the divider
through manystream.scan with exact images and exact partial
bounds; the root through a one-stream scan written here, since its
image is irrational: the box [x0, x1] maps to
[sqrt(s + x0), sqrt(s + x1)], only boxes with
4 Lam_law^2 (s + x0) < 1 can kill (the excess region, rational),
and the zone test asks for an integer strictly inside
(sqrt(s + x0)/e + M- - 1, sqrt(s + x1)/e - M+), e = b^(L - n),
decided by an exact floor of sqrt(A)/e plus a rational, by
squaring. The same one-stream scan is run on the reciprocal as
its control against manystream.scan.
THE LAW'S BAND over six families -- sqrt(s + y), 1/(s + y),
1/(s + y)^2, x/(s + y), exp(-h y), x exp(-h y) -- at the 20 cells
and h = j/(b - 1), j in {1, 2, 4, 8, 32, 128, 1024, 8192}; the
exponentials' ceilings in floating point, a pair whose argument sits
within 1e-9 of an integer printed as a near-tie and not read.
THE DIVIDER'S DENSE SWEEP: c* at h = j/(b - 1), j = 1..4096, and
h = q/10, q = 1..1000, at the 20 cells, by the closed forms.
THE LITERATURE'S TABLE: Ercegovac and Lang, Digital Arithmetic (2003),
Table 9.1, read in the authors' own chapter 9 slides, gives the
most-significant-first delays addition 2 (r = 2) and 1 (r >= 4),
multiplication 3 (r = 2) and 2 (r = 4), division 4, square root 4,
max/min 0; the slide deriving the multiplier's delay lists
(r, redundancy factor a/(r - 1), delta) = (2, 1, 3), (4, 1, 2),
(4, 2/3, 3). Their delta counts operand digits past the result
digit with operands and result in one format, so it is c at o = 0.
The product x y's c* at the three digit sets (2,1,1), (4,3,3) and
(4,2,2) is set beside it.
PREDICTIONS.
PA the controls, read first: the exact floor of sqrt against
60-digit decimals at 3000 random rationals, perfect squares
included, with no disagreement; the one-stream scan and
manystream.scan give the reciprocal the same verdict and the same
kill depth at every pair at L* - 1 and L*.
PB the closed forms equal the loops at all 360 (unit, pair); the
divider's c* is in {0, 1, 2} at every pair; the root's c* is
negative at 65 of its 120 pairs (a transplant from the root's
earlier census), the reciprocal's negative at some pair at every
cell.
PC at every (unit, pair) the scan at L* finds no kill and the scan
at L* - 1 finds one, the negative lookaheads included; the
deepest output depth at most 9.
PD the band holds at every (family, cell, h).
PE the divider: c* >= 0 at every dense pair, >= 1 at every slack-1
cell, <= 2 at every h >= 1/(b - 1); 0 at some h at every cell of
slack 2 or more (a transplant: the earlier sweep's zeros lay on
exactly those ten cells).
PF exp(-h y)'s c* at j = 8192 exceeds its j = 1 value at every cell;
the root's and both reciprocals' sit below theirs at every cell.
PG the product's c* is 2, 1 and 2 at (2,1,1), (4,3,3) and (4,2,2),
one below the literature's delta at each.
KILLS, as what the script prints.
K1 a PA line off: the comparator or the one-stream scan is wrong;
nothing below is read.
K2 a closed form off a loop: (2) or (3) is wrong.
K3 a kill at L*, or a pair alive at L* - 1 whose scan ran out of
region: the first is (1)'s sufficiency, the second necessity
unread there.
K4 the band broken at any read pair: (3) is wrong.
K5 a divider c* below 0 anywhere, or below 1 at slack 1: (4) is
wrong.
K6 exp(-h y) not rising at some cell: (6) is wrong.
POSITIVE CONTROL: PA, read before any verdict line.
THE ESTIMATE, a second question put after the record run: is the
literature's one extra multiplier digit (PG) the price of selecting the
output digit from an ESTIMATE of the residual rather than from the
residual itself?
THE LITERATURE'S CONDITION (Ercegovac and Lang, chapter 9 slides 43-44,
read at source, their radix, redundancy factor, estimate width and
delay written here b, Mh, m and L). On the symmetric set {-a, ...,
a} in radix b, with Mh = a/(b - 1), the multiplier's scaled
residual is bounded by om = Mh (1 - 2 Mh b^-L), and selection
constants exist when the residual is estimated to m fractional bits
iff T_m(om) >= (1 + 2^-m)/2, T_m the truncation to m fractional
bits. The slides tabulate (b, Mh, m, L) = (2, 1, 2, 3),
(4, 1, 2, 2) and (4, 2/3, 3, 3).
THE ARGUMENT (written before the engine). T_m(om) <= om and
(1 + 2^-m)/2 > 1/2, so some m works iff om > 1/2 STRICTLY: om > 1/2
gives T_m(om) > om - 2^-m, which clears (1 + 2^-m)/2 once
3 * 2^-m / 2 < om - 1/2. And om >= 1/2 unpacks to
b^L (2a + 1 - b)(b - 1) >= 4a^2, the product's law on the
symmetric set with lead 0. So the least delay the
condition reaches at any finite m is L* where the law holds
strictly at L*, and L* + 1 where it holds with EQUALITY. At
a = b - 1 equality reads b^L = 4: radix 2 at L = 2 and radix 4 at
L = 1, two of the three published cells.
PREDICTIONS.
PH1 (the control) the three published rows satisfy the condition,
and each fails at L - 1 at its own m.
PH2 the least L over m <= 64 is L* + 1 at (2,1,1) and (4,3,3),
where the law holds with equality at L*, and L* at (4,2,2),
first reached at m = 4.
PH3 over every symmetric cell of slack >= 1 at radices 2..64, the
least L over m <= 64 is L* where the law is strict and
L* + 1 where it is an equality; the equality cells are the two
of PH2 (a guess, not derived for a < r - 1). The formula's L*
equals manystream.law_of's c* at the symmetric cells of radix
<= 5.
KILLS, as prints.
K7 a PH1 row off: the condition is misread from the slides; every
verdict below is void (they still print).
K8 a cell whose least L is neither the argument's value nor
unreached at m = 64: the argument is wrong.
FINDINGS (entered after the run).
- PH1 holds: the three published rows satisfy the condition at
their L and fail at L - 1. K7 did not fire.
- PH2 holds: om at L* is exactly 1/2 at (2,1,1) and (4,3,3), where
the least L over every m is 3 and 2, the published values;
at (4,2,2) om = 11/18 and L = 2 is reached at m = 4, one
below the published 3 at m = 3.
- PH3 holds (the loop reads a <= b - 1, narrower than its design):
over the symmetric cells with a <= b - 1 of radices
2..64 the least L is L* at every strict cell and L* + 1 at the
equality cells, which are the two of PH2 and no others; the strict cells
reach L* by m = 8 (m = 2 at 460 cells, 3 at 273, 4 at 138, 5 at 72,
6 at 44, 7 at 33, 8 at 2). K8 did not fire.
So the extra digit is the estimate's only where the law is an
equality; elsewhere a wider estimate reaches the floor.
FINDINGS (entered after the run; every number is a print).
- PA holds: the exact floor of sqrt agrees with 60-digit decimals
at 3000 of 3000 rationals, and the one-stream scan with
manystream.scan on the reciprocal at 240 of 240 (pair, L). K1 did
not fire.
- PB holds: the closed forms equal the loops at 360 of 360. The
root's c* runs from -7 to 0 and is negative at 65 of 120 pairs,
as the transplant said; the reciprocal's from -6 to 2, negative at
40 pairs over all 20 cells; the divider's c* is 0 at 5 pairs, 1 at
75 and 2 at 40.
- PC holds: all 360 (unit, pair) are alive at L* and dead at
L* - 1, the 105 negative lookaheads included, every survival at
L* certified by the excess region running out. The kills sit at
output depth 1 to 9 for the root (the 9 at one pair), 1 to 8 for
the reciprocal and 1 to 8 for the divider.
- PD holds: the band at all 960 (family, cell, h), no near-tie. The
greatest deviations are 0.999 (the root), 0.876 and 0.893 (the
reciprocals), 1.000 printed for the divider (a float rounding of a
deviation under 1, which the check reads as under 1), 0.924
(exp(-h y)) and 0.952 (x exp(-h y)).
- PE holds: over 98620 dense pairs the divider's c* runs 0 to 4, at
most 2 at every h >= 1/(b - 1), at least 1 at every slack-1 cell,
and 0 somewhere at exactly the ten cells of slack 2 or more.
- PF holds: from j = 1 to j = 8192, exp(-h y)'s c* rises by 5 to 13
at every cell; the root's falls by 4 to 12, the reciprocal's and
1/(s + y)^2's by 5 to 13.
- PG holds: the product's c* is 2, 1 and 2 at (2,1,1), (4,3,3) and
(4,2,2), each dead at L* - 1, against the literature's delta 3, 2 and
3: the published multiplier sits one digit above the floor at all
three digit sets.
- PB's loops share each unit's Lam and range with its closed form,
so PB checks the ceilings and the divider's lead; PC reads Lam
against the scans.
- The reciprocal's PC scans ran through the one-stream scan, not
manystream.scan as the design says. PA shows that the two agree
at exactly those 240 (pair, L).
- One change after the first run, no verdict touched. The divider's
scan at L* - 1 as frozen (manystream.scan over both streams) ran
out of budget at (2,1,1), h = 128, where the excess is 0.008 of a
zone. Given more budget it found the kill at output depth 7 but
passed the 512 MB ceiling. Below the law the divider is now
searched on the edge |x| = Mh alone (edge_scan): a search for a
witness, which can only find a kill. Survival at L* is still read
by the full scan, whose region is empty at the root box.
RUN RECORD: 21/21 checks, 2.1 s, peak commit 14.4 MB under a memory guard;
after an audit (survival at L* read only where the excess region ran
out, PG's kills asserted), 21/21, 2.0 s, 13.8 MB.
"""
import math
import random
import time
from decimal import Decimal, getcontext
from fractions import Fraction as Fr
import manystream as ms
CHECKS = []
NMAX = 16
BUDGET = 300_000
JS = (1, 2, 4, 8, 32, 128)
JS_BAND = (1, 2, 4, 8, 32, 128, 1024, 8192)
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
tail = f" ({detail})" if detail else ""
print(f" {'PASS' if ok else 'FAIL'} {name}{tail}")
def least_pow(b, X):
"""The least integer k with b^k >= X, X a positive rational."""
k = math.log(X.numerator) - math.log(X.denominator)
k = math.floor(k / math.log(b)) - 2
while Fr(b) ** k < X:
k += 1
while Fr(b) ** (k - 1) >= X:
k -= 1
return k
# ---- exact values: a rational, or the square root of one ----
def sqrt_ge(A, y):
"""sqrt(A) >= y, exactly."""
return y <= 0 or y * y <= A
def vfloor(v, e, c):
"""floor(v/e + c), v = ('q', r) or ('r', A) meaning sqrt(A)."""
kind, x = v
if kind == "q":
return math.floor(x / e + c)
m = math.floor(math.sqrt(float(x)) / float(e) + float(c)) + 2
while not sqrt_ge(x, (m - c) * e):
m -= 1
return m
def vgt(v, e, y):
"""v/e > y, exactly."""
kind, x = v
if kind == "q":
return x / e > y
return y * e < 0 or (y * e) ** 2 < x
def one_kills(cell, lo, hi, n, L):
"""manystream's zone test on an image [lo, hi] of exact values."""
e = Fr(cell.b) ** (L - n)
return vgt(hi, e, vfloor(lo, e, cell.Mm) + cell.Mp)
# ---- the units: s = M- + h ----
class Root:
name = "sqrt(s + x)"
def __init__(self, cell, h):
self.c, self.h, self.s = cell, h, cell.Mm + h
def law_L(self):
X = Fr(self.c.I ** 2, 4 * self.c.rho ** 2) / self.h
return -((-least_pow(self.c.b, X)) // 2)
def lead(self):
X = (self.h + self.c.w) / self.c.Mp ** 2
return -((-least_pow(self.c.b, X)) // 2)
def loops(self):
c, L = self.c, -60
while 4 * Fr(c.b) ** (2 * L) * c.rho ** 2 * self.h < c.I ** 2:
L += 1
o = -60
while c.Mp ** 2 * Fr(c.b) ** (2 * o) < self.h + c.w:
o += 1
return L, o
def image(self, x0, x1):
return ("r", self.s + x0), ("r", self.s + x1)
def steep(self, x0, x1, law):
return 4 * law * law * (self.s + x0) < 1
def gamma(self):
h, w = float(self.h), float(self.c.w)
return float(self.c.Mp) / (2 * math.sqrt(h * (h + w)))
class Recip:
name = "1/(s + x)"
def __init__(self, cell, h):
self.c, self.h, self.s = cell, h, cell.Mm + h
self.lam = 1 / (h * h)
self.rng = (1 / (h + cell.w), 1 / h)
def law_L(self):
return least_pow(self.c.b, self.lam * self.c.I / self.c.rho)
def lead(self):
return least_pow(self.c.b, self.rng[1] / self.c.Mp)
def loops(self):
return self.c.law_L(self.lam), self.c.lead(self.rng)
def image(self, x0, x1):
return ("q", 1 / (self.s + x1)), ("q", 1 / (self.s + x0))
def steep(self, x0, x1, law):
return (self.s + x0) ** 2 * law < 1
def gamma(self):
return float(self.c.Mp / self.h)
def mapping(self):
s = self.s
return (self.name, 1,
lambda B: (1 / (s + B[0][1]), 1 / (s + B[0][0])),
lambda B: [(-1 / (s + B[0][0]) ** 2,
-1 / (s + B[0][1]) ** 2)],
None)
class Recip2:
name = "1/(s + y)^2"
def __init__(self, cell, h):
self.c, self.h = cell, h
def law_L(self):
return least_pow(self.c.b, 2 / self.h ** 3
* self.c.I / self.c.rho)
def lead(self):
return least_pow(self.c.b, 1 / (self.h ** 2 * self.c.Mp))
def gamma(self):
return float(2 * self.c.Mp / self.h)
class Divider:
name = "x/(s + y)"
def __init__(self, cell, h):
self.c, self.h, self.s = cell, h, cell.Mm + h
self.lam = (h + cell.Mh) / (h * h)
self.rng = (-cell.Mm / h, cell.Mp / h)
def law_L(self):
return least_pow(self.c.b, self.lam * self.c.I / self.c.rho)
def lead(self):
return least_pow(self.c.b, 1 / self.h)
def loops(self):
return self.c.law_L(self.lam), self.c.lead(self.rng)
def gamma(self):
return float(1 + self.c.Mh / self.h)
def mapping(self):
s = self.s
def image(B):
(x0, x1), (y0, y1) = B
v = [x / (s + y) for x in (x0, x1) for y in (y0, y1)]
return (min(v), max(v))
def partials(B):
(x0, x1), (y0, y1) = B
X = max(abs(x0), abs(x1))
return [(1 / (s + y1), 1 / (s + y0)),
(-X / (s + y0) ** 2, X / (s + y0) ** 2)]
return (self.name, 2, image, partials, None)
def one_scan(cell, unit, L, o):
"""The one-stream scan: least output depth of a kill at L, or None,
with the boxes read and whether the excess region ran out."""
c = L - o
law = cell.lam_law(L)
D = range(-cell.am, cell.ap + 1)
level, read = [0], 0
for n in range(NMAX + 1):
kept = []
for u in level:
x0, x1 = cell.box(u, n)
read += 1
if not unit.steep(x0, x1, law):
continue
lo, hi = unit.image(x0, x1)
if n - c >= 1 and one_kills(cell, lo, hi, n, L):
return n - c, read, False
kept.append(u)
if not kept:
return None, read, True
if read > BUDGET:
return None, read, False
level = [cell.b * u + e for u in kept for e in D]
return None, read, False
def many_scan(cell, unit, L, o):
return ms.scan(cell, unit.mapping(), L, o)
def edge_scan(cell, unit, L, o):
"""A witness search for the divider below its law: the x prefix
held on the edge |x| = Mh, only the y prefixes branching. A kill
found here is a kill; a survival here reads nothing."""
_, _, image, partials, _ = unit.mapping()
c, law = L - o, cell.lam_law(L)
top = cell.Mp >= cell.Mm
level, read = [0], 0
for n in range(NMAX + 1):
ux = ms.repunit(cell.b, n) * (cell.ap if top else -cell.am)
kept = []
for uy in level:
B = [cell.box(ux, n), cell.box(uy, n)]
read += 1
if ms.rate(partials(B)) <= law:
continue
if n - c >= 1 and cell.kills(image(B), n, L):
return n - c, read
kept.append(uy)
if not kept or read > BUDGET:
return None, read
level = [cell.b * u + e for u in kept
for e in range(-cell.am, cell.ap + 1)]
return None, read
# ---- the sections ----
def pairs():
for (b, am, ap) in ms.census():
cell = ms.Cell(b, am, ap)
for j in JS:
yield cell, Fr(j, b - 1)
def section_pa():
print("PA the controls")
getcontext().prec = 60
rng = random.Random(7)
bad = 0
for i in range(3000):
if i % 3 == 0:
q = Fr(rng.randint(1, 400), rng.randint(1, 60))
A = q * q
else:
A = Fr(rng.randint(1, 10 ** 6), rng.randint(1, 10 ** 4))
e = Fr(rng.randint(1, 50), rng.randint(1, 50))
c = Fr(rng.randint(-30, 30), rng.randint(1, 9))
exact = vfloor(("r", A), e, c)
z = (Decimal(A.numerator) / Decimal(A.denominator)).sqrt()
z = z / (Decimal(e.numerator) / Decimal(e.denominator))
z += Decimal(c.numerator) / Decimal(c.denominator)
if exact != int(z.to_integral_value(rounding="ROUND_FLOOR")):
bad += 1
check("PA the exact floor of sqrt against 60-digit decimals",
bad == 0, f"3000 rationals, {bad} off")
off = tot = 0
for cell, h in pairs():
u = Recip(cell, h)
Ls, o = u.law_L(), u.lead()
for L in (Ls - 1, Ls):
k1, _, _ = one_scan(cell, u, L, o)
k2, _ = many_scan(cell, u, L, o)
tot += 1
off += k1 != k2
check("PA the one-stream scan equals manystream.scan on 1/(s + x)",
off == 0, f"{tot} (pair, L), {off} off")
def section_pb_pc():
print("PB, PC the three units at the 120 pairs")
t0 = time.time()
loops_off, neg, cstar, negcells = 0, {}, {}, {}
alive_bad, dead_bad, depth, reads, least = [], [], {}, 0, {}
for U in (Root, Recip, Divider):
neg[U.name], cstar[U.name], depth[U.name] = 0, {}, {}
negcells[U.name] = set()
for cell, h in pairs():
u = U(cell, h)
Ls, o = u.law_L(), u.lead()
loops_off += (Ls, o) != u.loops()
cs = Ls - o
cstar[U.name][cs] = cstar[U.name].get(cs, 0) + 1
at = (cs, (cell.b, cell.am, cell.ap, str(h)))
least[U.name] = min(least.get(U.name, at), at)
if cs < 0:
neg[U.name] += 1
negcells[U.name].add((cell.b, cell.am, cell.ap))
for L in (Ls, Ls - 1):
if U is Divider and L < Ls:
k, r = edge_scan(cell, u, L, o)
ran_out = False
elif U is Divider:
k, r = many_scan(cell, u, L, o)
ran_out = (k is None and r <= ms.BUDGET
and ms.scan.reached < ms.NMAX)
else:
k, r, ran_out = one_scan(cell, u, L, o)
reads += r
tag = (U.name, cell.b, cell.am, cell.ap, str(h), L)
# alive at L* only where the excess region ran out
if L == Ls and (k is not None or not ran_out):
alive_bad.append(tag + (ran_out,))
if L == Ls - 1:
if k is None:
dead_bad.append(tag + (ran_out,))
else:
depth[U.name][k] = depth[U.name].get(k, 0) + 1
print(f" {U.name}: c* {dict(sorted(cstar[U.name].items()))}, "
f"negative at {neg[U.name]} pairs over "
f"{len(negcells[U.name])} cells, least {least[U.name][0]} at "
f"{least[U.name][1]}; kills at L* - 1 by "
f"output depth {dict(sorted(depth[U.name].items()))}")
check("PB the closed forms equal the loops", loops_off == 0,
f"360 (unit, pair), {loops_off} off")
check("PB the divider's c* in {0, 1, 2}",
set(cstar["x/(s + y)"]) <= {0, 1, 2})
check("PB the root negative at 65 of 120 pairs",
neg["sqrt(s + x)"] == 65, f"read {neg['sqrt(s + x)']}")
check("PB the reciprocal negative at every cell",
len(negcells["1/(s + x)"]) == 20,
f"read {len(negcells['1/(s + x)'])} cells")
deepest = max(max(d) for d in depth.values() if d)
check("PC alive at L*, dead at L* - 1, at every (unit, pair)",
not alive_bad and not dead_bad,
f"killed at L*: {alive_bad}; alive at L* - 1: {dead_bad}")
check("PC the deepest kill at output depth <= 9", deepest <= 9,
f"deepest {deepest}")
print(f" {reads} boxes read, {time.time() - t0:.1f} s")
def exp_c(cell, h, times_x):
"""c* and the band deviation of exp(-h y) (times x if asked), in
floating point; None when a ceiling's argument is a near-tie."""
lb = math.log(cell.b)
Mm, Mp, Mh = float(cell.Mm), float(cell.Mp), float(cell.Mh)
hf = float(h)
ratio = math.log(cell.I / cell.rho)
if times_x:
lam = math.log(1 + Mh * hf) + hf * Mm
sig = hf * Mm
gam = 1 + Mh * hf
else:
lam = math.log(hf) + hf * Mm
sig = hf * Mm - math.log(Mp)
gam = hf * Mp
X, Y = (lam + ratio) / lb, sig / lb
if min(abs(X - round(X)), abs(Y - round(Y))) < 1e-9:
return None
cs = math.ceil(X) - math.ceil(Y)
return cs, cs - math.log(gam * cell.I / cell.rho) / lb
def section_pd_pf():
print("PD, PF the relative-gradient law over six families")
worst, broken, ties, read = {}, [], 0, 0
first, last = {}, {}
fams = [Root, Recip, Recip2, Divider, "exp(-h y)", "x exp(-h y)"]
for (b, am, ap) in ms.census():
cell = ms.Cell(b, am, ap)
for j in JS_BAND:
h = Fr(j, b - 1)
for F in fams:
if isinstance(F, str):
got = exp_c(cell, h, F.startswith("x"))
if got is None:
ties += 1
continue
cs, dev = got
name = F
else:
u = F(cell, h)
name = u.name
cs = u.law_L() - u.lead()
dev = cs - math.log(u.gamma() * cell.I
/ cell.rho) / math.log(b)
worst[name] = max(worst.get(name, 0), abs(dev))
read += 1
if not abs(dev) < 1:
broken.append((name, b, am, ap, j))
if j == JS_BAND[0]:
first[name, b, am, ap] = cs
if j == JS_BAND[-1]:
last[name, b, am, ap] = cs
print(" greatest |deviation|: " + ", ".join(
f"{k} {v:.3f}" for k, v in worst.items()) + f"; near-ties {ties}; "
f"{read} (family, cell, h) read")
check("PD the band holds at every read pair", not broken,
f"broken at {broken}")
for name, sign in (("exp(-h y)", 1), ("sqrt(s + x)", -1),
("1/(s + x)", -1), ("1/(s + y)^2", -1)):
moves = [sign * (last[k] - first[k]) for k in first
if k[0] == name and k in last]
print(f" {name}: c* at j = 8192 less c* at j = 1 runs "
f"{min(sign * m for m in moves)} to "
f"{max(sign * m for m in moves)}")
check(f"PF {name} {'rises' if sign > 0 else 'falls'} "
f"at every cell", len(moves) == 20 and min(moves) > 0)
def section_pe():
print("PE the divider's dense sweep")
t0 = time.time()
seen, low, high, zero_cells = 0, [], [], set()
lo_all, hi_all = 99, -99
for (b, am, ap) in ms.census():
cell = ms.Cell(b, am, ap)
hs = {Fr(j, b - 1) for j in range(1, 4097)}
hs |= {Fr(q, 10) for q in range(1, 1001)}
for h in hs:
u = Divider(cell, h)
cs = u.law_L() - u.lead()
seen += 1
lo_all, hi_all = min(lo_all, cs), max(hi_all, cs)
if cs < 0 or (cell.rho == 1 and cs < 1):
low.append((b, am, ap, str(h), cs))
if h >= Fr(1, b - 1) and cs > 2:
high.append((b, am, ap, str(h), cs))
if cs == 0:
zero_cells.add((b, am, ap))
slack2 = {(b, am, ap) for (b, am, ap) in ms.census()
if am + ap + 1 - b >= 2}
print(f" {seen} pairs, c* from {lo_all} to {hi_all}; zeros at "
f"{len(zero_cells)} cells; {time.time() - t0:.1f} s")
check("PE c* >= 0 everywhere, >= 1 at slack 1", not low,
f"off at {low[:5]}")
check("PE c* <= 2 at every h >= 1/(b - 1)", not high,
f"off at {high[:5]}")
check("PE a zero at every cell of slack 2 or more, none at slack 1",
zero_cells == slack2,
f"{len(zero_cells)} cells, {len(slack2)} of slack >= 2")
def section_pg():
print("PG the product against the literature's multiplier delays")
xy = ms.MAPS[2]
got, kills = [], []
for (b, am, ap), delta in (((2, 1, 1), 3), ((4, 3, 3), 2),
((4, 2, 2), 3)):
cell = ms.Cell(b, am, ap)
lam, o, Ls = ms.law_of(cell, xy)
k, _ = ms.scan(cell, xy, Ls - 1, o)
got.append(Ls - o)
kills.append(k)
print(f" ({b}, {am}, {ap}): lead {o}, L* {Ls}, c* {Ls - o}, "
f"dead at L* - 1 at output depth {k}; the literature's delta "
f"{delta}")
check("PG c* = 2, 1, 2, each dead at L* - 1",
got == [2, 1, 2] and None not in kills, f"read {got}")
def trunc(x, t):
"""x truncated to t fractional bits."""
return Fr(math.floor(x * 2 ** t), 2 ** t)
def el_om(r, a, delta):
p = Fr(a, r - 1)
return p * (1 - 2 * p / Fr(r) ** delta)
def el_ok(r, a, t, delta):
return trunc(el_om(r, a, delta), t) >= Fr(1 + Fr(1, 2 ** t), 2)
def el_least(r, a, tmax=64, dmax=40):
"""(least delta over t <= tmax, the least t reaching it)."""
for delta in range(0, dmax):
for t in range(1, tmax + 1):
if el_ok(r, a, t, delta):
return delta, t
return None, None
def xy_law(r, a):
"""(L*, equality at L*) of r^L (2a + 1 - r)(r - 1) >= 4a^2."""
L = 0
while r ** L * (2 * a + 1 - r) * (r - 1) < 4 * a * a:
L += 1
return L, r ** L * (2 * a + 1 - r) * (r - 1) == 4 * a * a
def section_ph():
print("PH the literature's selection condition at an exact estimate")
rows = ((2, 1, 2, 3), (4, 3, 2, 2), (4, 2, 3, 3))
bad = [(r, a, t, d) for r, a, t, d in rows
if not el_ok(r, a, t, d) or el_ok(r, a, t, d - 1)]
check("PH1 the published rows hold at L and fail at L - 1",
not bad, f"off at {bad}")
for r, a, _, d in rows:
L, eq = xy_law(r, a)
dl, tl = el_least(r, a)
print(f" ({r}, {a}, {a}): om at L* {el_om(r, a, L)}, law "
f"L* {L}{' with equality' if eq else ''}; least L "
f"{dl} first at m = {tl}; published {d}")
got = [el_least(r, a) for r, a, _, _ in rows]
check("PH2 least L 3, 2, 2, the last first at m = 4",
[g[0] for g in got] == [3, 2, 2] and got[2][1] == 4,
f"read {got}")
off, unreached, eqs, tneed = [], [], [], {}
for r in range(2, 65):
for a in range((r + 1) // 2, r):
L, eq = xy_law(r, a)
dl, tl = el_least(r, a)
if dl is None:
unreached.append((r, a))
elif dl != L + (1 if eq else 0):
off.append((r, a, L, eq, dl))
if eq:
eqs.append((r, a, L))
if dl == L:
tneed[tl] = tneed.get(tl, 0) + 1
print(f" radices 2..64: equality cells (b, a, L) {eqs}; least m reaching "
f"L* at the strict cells {dict(sorted(tneed.items()))}; "
f"unreached {unreached}")
check("PH3 least L = L* strict, L* + 1 at equality", not off,
f"off at {off[:6]}")
lawoff = []
for (b, am, ap) in ms.census():
if am == ap:
cell = ms.Cell(b, am, ap)
lam, o, Ls = ms.law_of(cell, ms.MAPS[2])
if Ls - o != xy_law(b, am)[0]:
lawoff.append((b, am))
check("PH3 the formula's L* equals manystream's c* at radix <= 5",
not lawoff, f"off at {lawoff}")
def main():
t0 = time.time()
section_pa()
if not all(CHECKS):
print("\nCONTROL FAILED: nothing below is read.")
raise SystemExit(1)
section_pb_pc()
section_pd_pf()
section_pe()
section_pg()
section_ph()
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed, "
f"{time.time() - t0:.1f} s")
raise SystemExit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()