# Python program for OEIS A290154 # Michael S. Branicky, Jun 20 2021 # (Python) from sympy import factorint, prevprime, primerange, prod def aupto(limit, bfile=False): adict, pN = dict(), prevprime(limit+1) pi = {p: i for i, p in enumerate(primerange(1, pN+1), start=1)} smooth = {i: 0 for i in pi.values()} watching = smooth[1] = 1 # 1 is prime(n) smooth for all n for n in range(2, limit+1): f = factorint(n, limit=pN) nt = prod(p**f[p] for p in f if p <= pN) if nt == n: bigpi = pi[max(f)] if bigpi <= watching: smooth[watching] += 1 else: smooth[bigpi] += 1 if 2*smooth[watching] == n: adict[watching] = n if bfile: print(watching, n) with open('b290154.txt', 'a') as bfile: bfile.write(f"{watching} {n}\n") if watching == 10000: break watching += 1 smooth[watching] += smooth[watching-1] return sorted(adict.values()) print(aupto(12000)) # ~~~~ aupto(3*10**8, bfile=True) # ~~~~