"""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())