dvalue.py
Checks: HEADROOM: The blind-density theorem theorem; rule.
Cited on: HEADROOM.
To run it: save dvalue.py, with growth.py, headroom.py beside it; needs numpy.
"""dvalue.py -- the value of the blind density at lambda = 2: how do the
truncations delta_T reach their limit?
QUESTION. At lambda = 2 the class of W(2) = 24 is the primes q = 2n + 1
with n coprime to 6 and no divisor 1 < d < n of n with 2d + 1 prime
(q = 11 mod 12 is v_2(q - 1) <= 1 with 3 not dividing q - 1, the door
3; p - 1 | 2n forces
(p - 1)/2 | n, which is coprime to 6, and a divisor d = 1 mod 6 makes
2d + 1 divisible by 3). Under the prime density n's valuations are
independent, P(l^a | n) = 1/(l^(a-1) (l - 1)) for l >= 5, and
delta_T = (1/4) P(no such d <= (T - 1)/2),
nonincreasing in T, with a positive limit delta. Can a truncation read
delta to three digits?
DESIGN.
(1) THE UPPER BOUND (proved). The T-event reads the coordinates
l <= (T - 1)/2; "no prime l in ((T - 1)/2, X] with 2l + 1 prime divides
n" reads disjoint ones, and both contain the T'-event, T' = 2X + 1, so
delta <= delta_T' <= delta_T * prod (1 - 1/(l - 1)) over those l, the
Sophie Germain primes past (T - 1)/2, a product that converges (their sum of 1/l is finite by
the upper-bound sieve). Printed at T = 47 with delta_47 exact.
(2) THE MONTE CARLO. N samples of n's valuations at every prime 5 <= l
<= L = (TMAX - 1)/2: below 2000 an exact Bernoulli draw per prime;
above, a binomial count per prime and uniform sample indices, a sample
drawn twice for one l kept once (relative bias about count/2N, under
3 10^-4). The exponent is geometric, P(a >= k + 1 | a >= 1) = l^-k. Per
sample the least divisor d > 1, d <= L, with 2d + 1 prime is found, so
every delta_T, T <= TMAX, is read off one sample set, nested, and a
decade's increment is a count.
(3) THE PRIME SHARE. Over the primes q <= X = TMAX, mark q = 1 + k(p - 1),
k >= 2, for p = 5, 7 and every prime p = 11 mod 12 (another p is caught
by 5 or 7: 4 | p - 1 or 6 | p - 1); S(x) is the share of the primes
<= x left unmarked, q >= 11, and S_T(x) the same with p <= T only.
PREDICTIONS (frozen before the code).
R1 (property) delta_T nonincreasing; delta_47 = 0.16108 exactly.
R2 (pattern, 70%) the decade increments I_j = delta_(10^j) -
delta_(10^(j+1)) shrink slowly: I_5/I_6 and I_6/I_7 below 1.5. A
tail (log T)^(-c), c <= 1, gives ratios ((j+1)/j)^(1+c) <= 1.36 at
j = 6; a tail T^(-a), a >= 0.2, gives 10^a >= 1.58.
R3 (guess, 55%) delta_(10^8) < 0.150.
R4 (TRANSPLANT, re-derive on the print) S(x) sits within 0.003 of
delta_x at 10^8.
THE KILL (of "a truncation reads the value"), printed: "LOGARITHMIC"
when both ratios I_5/I_6, I_6/I_7 < 1.5, "SETTLING" when both >= 1.5,
else "MIXED"; a ratio is read only when each count behind it is >= 400
(5% error), else "UNREAD".
CONTROLS, read first. M1 the Monte Carlo at T = 11, 23, 47 within 4
sigma of the exact delta_T (headroom.delta_exact). M2 planted: draw
P(l | n) = 1/l (the integers, not shifted primes); M1 must fail at 47.
M3 the prime share with 5 <= p <= 47 only, over the primes in (47,
10^7], prints headroom's 0.16100.
FINDINGS (from the printed run; the controls passed first). M1: the
Monte Carlo reads delta_11, delta_23, delta_47 as 0.18759, 0.16886,
0.16113 against the exact 0.18750, 0.16875, 0.16108 (z = +0.79, +0.95,
+0.42); M2: the planted 1/l reads 0.17402, z = +112.5; M3: 0.16100.
The truncations, each +- 0.00012:
T 47 10^2 10^3 10^4 10^5 10^6
delta_T 0.16113 0.15156 0.13283 0.12394 0.11792 0.11350
T 10^7 10^8
delta_T 0.11015 0.10750
The decade increments I_2..I_7 are 0.01873, 0.00889, 0.00602, 0.00442,
0.00335, 0.00265 (74912, 35567, 24061, 17677, 13401, 10617 samples),
their ratios 2.106, 1.478, 1.361, 1.319, 1.262: R2 holds, and THE KILL
prints LOGARITHMIC: no truncation reads the value. log I_j against
log j has slope -1.530: a pure tail (log T)^(-c) would have c near
0.53, ((j+1)/j)^1.53 giving 1.322 and 1.266 at j = 5, 6, but a shifted
tail K/(ln T + 3) gives 1.465, 1.377, 1.317, 1.274 at j = 3..6 as well
(by hand), so c is a model's, not the run's. The kill could not certify
three digits in any case: I_7 alone is 0.00265. R3 holds
(0.10750 < 0.150). The upper bound (1): delta <= 0.11969, the Sophie
Germain primes to 5 10^7 taken. The prime share S(x) at 10^4..10^8 is
0.12937, 0.11812, 0.11172, 0.10747, 0.10445, below delta_x by 0.00178,
0.00269, 0.00305 at 10^6, 10^7, 10^8: R4 UNDECIDED at its bar, 0.00305
+- 0.00012 against 0.003, while the gap's growth, 0.00126 +- 0.00004
over the 24018 nested samples behind it, is read: the transplant does
not hold as x grows, unexplained (a candidate: an integer q - 1 <= x
pairs each divisor d with (q - 1)/d, which independent valuations do
not). S_T(10^8) with
p <= T, T = 10^2..10^8: 0.15146, 0.13273, 0.12308, 0.11553, 0.10941,
0.10503, 0.10445, within 0.0001 of delta_T at 10^2 and 10^3 and
falling away above. So the value of delta is not read by any
truncation run here: it lies in (0, 0.11969], proved, and the
truncations are still falling at 10^8 with decade ratios under 1.5 from
I_3/I_4 on,
slower than any power tail T^(-a) with a >= 0.18.
RUN RECORD. Rehearsal --tmax 10^6 --n 10^5: 3/4 (M3 needs 10^7), 0.4
s, peak 59.0 MB. A first full run was refused at the 512 MB ceiling
(the per-sample loop's lists built whole); chunked by 10^5 samples, the
full run (TMAX = 10^8, N = 10^6, seed 1) 4/4, 4.1 s, peak 343.7 MB
under a 512 MB commit ceiling, the gap line added after the first record and every
figure above unchanged by it.
python dvalue.py [--tmax T] [--n N] [--seed S] (10^8, 10^6, 1)
"""
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np # noqa: E402
import headroom # noqa: E402
RESULTS = []
def check(name, ok, detail=""):
RESULTS.append(bool(ok))
print(f" [{'PASS' if ok else 'FAIL'}] {name}" + (f" -- {detail}" if detail
else ""))
def sieve(n):
s = np.ones(n + 1, dtype=bool)
s[:2] = False
s[4::2] = False
for i in range(3, int(n ** 0.5) + 1, 2):
if s[i]:
s[i * i::2 * i] = False
return s
def unmarked(isp, bad, lo, hi):
"""Primes in [lo, hi] left unmarked, and all primes there, counted in
chunks so no array of the sieve's size is made."""
k = t = 0
for a in range(lo, hi + 1, 10 ** 7):
b = min(hi + 1, a + 10 ** 7)
k += int(np.count_nonzero(isp[a:b] & ~bad[a:b]))
t += int(np.count_nonzero(isp[a:b]))
return k, t
def draw(rng, ls, N, pinv):
"""Sample indices and exponents for each prime in ls; pinv(l) is the
chance that l divides n. Returns (idx, l, a) arrays, deduplicated."""
I, Ls, A = [], [], []
small = ls[ls < 2000]
for l in small:
hit = np.nonzero(rng.random(N) < pinv(int(l)))[0]
I.append(hit)
Ls.append(np.full(hit.size, l, dtype=np.int64))
big = ls[ls >= 2000]
if big.size:
cnt = rng.binomial(N, np.array([pinv(int(l)) for l in big]))
I.append(rng.integers(0, N, size=int(cnt.sum())))
Ls.append(np.repeat(big, cnt))
idx = np.concatenate(I).astype(np.int64)
lv = np.concatenate(Ls)
key = np.unique(idx * (int(lv.max()) + 1) + lv)
idx, lv = key // (int(lv.max()) + 1), key % (int(lv.max()) + 1)
a = rng.geometric(1.0 - 1.0 / lv)
return idx, lv, a
def least_d(idx, lv, a, N, L, isp):
"""Per sample, the least divisor d > 1, d <= L, of n with 2d + 1
prime; L + 1 when there is none."""
order = np.argsort(idx, kind="stable")
idx, lv, a = idx[order], lv[order], a[order]
out = np.full(N, L + 1, dtype=np.int64)
cuts = np.searchsorted(idx, np.arange(0, N + 10 ** 5, 10 ** 5))
for c0, c1 in zip(cuts[:-1].tolist(), cuts[1:].tolist()):
ci, cl, ca = idx[c0:c1], lv[c0:c1], a[c0:c1]
if not ci.size:
continue
starts = np.flatnonzero(np.r_[True, ci[1:] != ci[:-1]])
ends = np.r_[starts[1:], ci.size]
idx_l, lv_l, a_l = ci.tolist(), cl.tolist(), ca.tolist()
for s, e in zip(starts.tolist(), ends.tolist()):
divs = [1]
for j in range(s, e):
l, ex = lv_l[j], a_l[j]
new = []
for d in divs:
m = d
for _ in range(ex):
m *= l
if m > L:
break
new.append(m)
divs += new
best = L + 1
for d in divs:
if 1 < d < best and isp[2 * d + 1]:
best = d
out[idx_l[s]] = best
return out
def main():
args = sys.argv[1:]
tmax = int(float(args[args.index("--tmax") + 1])) if "--tmax" in args \
else 10 ** 8
N = int(float(args[args.index("--n") + 1])) if "--n" in args else 10 ** 6
seed = int(args[args.index("--seed") + 1]) if "--seed" in args else 1
t0 = time.time()
isp = sieve(tmax + 1)
print(f"sieve to {tmax + 1}: {int(isp.sum())} primes "
f"({time.time() - t0:.1f} s)")
print("THE PRIME SHARE")
primes = np.flatnonzero(isp)
bad = np.zeros(tmax + 1, dtype=bool)
marks = primes[(primes == 5) | (primes == 7)
| ((primes % 12 == 11) & (primes >= 11))]
marks = marks[2 * (marks - 1) + 1 <= tmax]
decades = [10 ** j for j in range(2, 9) if 10 ** j <= tmax]
st = {}
m3 = None
for p in marks.tolist():
if p > 47 and m3 is None and tmax >= 10 ** 7:
k, t = unmarked(isp, bad, 48, 10 ** 7)
m3 = k / t
while decades and p > decades[0]:
k, t = unmarked(isp, bad, 11, tmax)
st[decades.pop(0)] = k / (t + 4)
bad[1 + 2 * (p - 1)::p - 1] = True
for T in decades:
k, t = unmarked(isp, bad, 11, tmax)
st[T] = k / (t + 4)
sx = {}
for x in (10 ** j for j in range(4, 9)):
if x <= tmax:
k, t = unmarked(isp, bad, 11, x)
sx[x] = k / (t + 4)
del bad
print(" S(x), every p: " + ", ".join(
f"10^{len(str(x)) - 1} {v:.5f}" for x, v in sx.items()))
print(" S_T(X), p <= T: " + ", ".join(
f"10^{len(str(T)) - 1} {v:.5f}" for T, v in st.items()))
if m3 is not None:
print(f" S_47 over the primes in (47, 10^7]: {m3:.5f}")
print("CONTROLS")
L = (tmax - 1) // 2
ls = primes[(primes >= 5) & (primes <= L)].astype(np.int64)
del primes
rng = np.random.default_rng(seed)
t1 = time.time()
idx, lv, a = draw(rng, ls, N, lambda l: 1.0 / (l - 1))
dmin = least_d(idx, lv, a, N, L, isp)
print(f" Monte Carlo: {N} samples, {idx.size} prime factors at "
f"5 <= l <= {L} ({time.time() - t1:.1f} s)")
del idx, lv, a
exact = {}
ok1 = True
for T in (11, 23, 47):
ps = [p for p in range(5, T + 1) if isp[p]]
exact[T] = headroom.delta_exact(ps)
P = float(np.mean(dmin > (T - 1) // 2))
sig = 0.25 * (P * (1 - P) / N) ** 0.5
z = (0.25 * P - float(exact[T])) / sig
ok1 &= abs(z) < 4
print(f" T = {T}: exact {float(exact[T]):.5f}, Monte Carlo "
f"{0.25 * P:.5f}, z = {z:+.2f}")
check("M1 the Monte Carlo within 4 sigma of exact delta_T, T = 11, 23, "
"47", ok1)
prng = np.random.default_rng(seed + 1)
pls = ls[ls <= 23]
pi, pl, pa = draw(prng, pls, N, lambda l: 1.0 / l)
pd = least_d(pi, pl, pa, N, 23, isp)
P = float(np.mean(pd > 23))
sig = 0.25 * (P * (1 - P) / N) ** 0.5
z2 = (0.25 * P - float(exact[47])) / sig
check("M2 planted P(l | n) = 1/l fails M1 at 47", abs(z2) >= 4,
f"{0.25 * P:.5f}, z = {z2:+.1f}")
check("M3 the prime share at p <= 47 over (47, 10^7] prints 0.16100",
m3 is not None and f"{m3:.5f}" == "0.16100",
f"{m3:.5f}" if m3 is not None else "tmax below 10^7")
check("R1 delta_47 exact is 0.16108", f"{float(exact[47]):.5f}" ==
"0.16108")
print("THE TRUNCATIONS")
Ts = [47] + [10 ** j for j in range(2, 9) if 10 ** j <= tmax]
dT = {T: 0.25 * float(np.mean(dmin > (T - 1) // 2)) for T in Ts}
for T in Ts:
P = dT[T] / 0.25
print(f" delta_{T:<10d} {dT[T]:.5f} +- "
f"{0.25 * (P * (1 - P) / N) ** 0.5:.5f} S_T(X) "
f"{st.get(T, float('nan')):.5f}")
print(f" never by l <= {L}: {0.25 * float(np.mean(dmin > L)):.5f}")
xs = [x for x in (10 ** 6, 10 ** 7, 10 ** 8) if x in sx]
if len(xs) == 3:
c = int(np.sum((dmin > (10 ** 6 - 1) // 2)
& (dmin <= (10 ** 8 - 1) // 2)))
g = [dT[x] - sx[x] for x in xs]
print(" delta_x - S(x): " + ", ".join(
f"10^{len(str(x)) - 1} {v:.5f}" for x, v in zip(xs, g))
+ f"; growth {g[2] - g[0]:.5f} +- "
f"{0.25 * (c / N * (1 - c / N) / N) ** 0.5:.5f} ({c} samples)")
inc = {}
for j in range(2, 8):
lo, hi = (10 ** j - 1) // 2, (10 ** (j + 1) - 1) // 2
if 10 ** (j + 1) > tmax:
break
c = int(np.sum((dmin > lo) & (dmin <= hi)))
inc[j] = c
print(f" I_{j}: {c} samples first fail in (10^{j}, 10^{j + 1}], "
f"{0.25 * c / N:.5f}")
ratios = {}
for j in sorted(inc):
if j + 1 in inc:
ratios[j] = (inc[j] / inc[j + 1] if min(inc[j], inc[j + 1])
>= 400 else None)
print(" ratios I_j/I_(j+1): " + ", ".join(
f"{j}: {r:.3f}" if r else f"{j}: UNREAD" for j, r in ratios.items()))
js = [j for j in sorted(inc) if inc[j] > 0]
if len(js) >= 3:
slope = np.polyfit(np.log(js), np.log([inc[j] for j in js]), 1)[0]
print(f" log I_j against log j: slope {slope:.3f} (a tail "
f"(log T)^-c has slope -(1 + c))")
print("THE UPPER BOUND")
sg = 1.0
for l in ls[ls > 23].tolist():
if isp[2 * l + 1]:
sg *= 1.0 - 1.0 / (l - 1)
print(f" delta <= delta_47 * prod over Sophie Germain l in (23, {L}] "
f"= {float(exact[47]) * sg:.5f}")
print("THE KILL")
r5, r6 = ratios.get(5), ratios.get(6)
if r5 is None or r6 is None:
print(" UNREAD")
elif r5 < 1.5 and r6 < 1.5:
print(" LOGARITHMIC: no truncation reads the value")
elif r5 >= 1.5 and r6 >= 1.5:
print(f" SETTLING: delta = {dT[Ts[-1]]:.4f} within about I_7")
else:
print(" MIXED")
print(f"{sum(RESULTS)}/{len(RESULTS)} checks passed, "
f"{time.time() - t0:.1f} s")
return 0 if all(RESULTS) else 1
if __name__ == "__main__":
sys.exit(main())