# Python program for OEIS A254859 # Michael S. Branicky, May 25 2022 # A254859 Numbers that are both a sum and a product of two or more consecutive primes. data = [15, 30, 77, 143, 210, 221, 323, 1001, 2310, 4199, 5767, 7429, 9797, 10403, 11021, 12317, 20711, 22499, 23707, 25591, 28891, 30030, 33263, 34571, 36863, 38021, 46189, 47053, 75067, 77837, 79523, 82861, 82919, 89951, 95477, 99221, 104927, 111547, 116939, 136891, 141367, 145157, 146969, 154433] import heapq from sympy import sieve def prime(n): return sieve[n] def prodgen(): p = prime(1) * prime(2) h = [(p, 1, 2)] nextcount = 3 while True: (v, s, l) = heapq.heappop(h) yield v if v >= p: p *= prime(nextcount) heapq.heappush(h, (p, 1, nextcount)) nextcount += 1 v //= prime(s); s += 1; l += 1; v *= prime(l) heapq.heappush(h, (v, s, l)) def sumgen(): p = prime(1) + prime(2) h = [(p, 1, 2)] nextcount = 3 oldv = -1 while True: (v, s, l) = heapq.heappop(h) if v != oldv: yield v oldv = v if v >= p: p += prime(nextcount) heapq.heappush(h, (p, 1, nextcount)) nextcount += 1 v -= prime(s); s += 1; l += 1; v += prime(l) heapq.heappush(h, (v, s, l)) def agen(): gp = prodgen() kp = next(gp) for ks in sumgen(): while kp < ks: kp = next(gp) if ks == kp: yield kp g = agen() ans = [next(g) for n in range(1, 1+len(data))] print(ans) print(data) assert ans == data print() from time import time time0 = time() for n, an in enumerate(agen(), 1): with open('b254859.txt', 'a') as bfile: bfile.write(f"{n} {an}\n") if n%100 == 0: print(n, an, time()-time0)