# Python program for OEIS AA059958 # after my translation of Daniel Suteu's PARI for 335645 # Michael S. Branicky, Feb 08 2023 # A059958 Smallest number m such that m*(m+1) has at least n distinct prime factors. 4 data =[1, 2, 5, 14, 65, 209, 714, 7314, 38570, 254540, 728364, 11243154, 58524465, 812646120, 5163068910, 58720148850, 555409903685, 4339149420605] from sympy import primorial, primerange, integer_nthroot, isprime, factorint def cond(n): m = integer_nthroot(n, 2)[0] return m*(m+1) == n def omega_cond(A, B, n): A = max(A, primorial(n)) def f(m, p, j): lst = [] for q in primerange(p, integer_nthroot(B//m, max(j, 1))[0]+2): v = m*q while v <= B: if j <= 1: if v >= A and cond(v): lst.append(v) elif v*(q+1) <= B: lst += f(v, q+1, j-1) v *= q return lst return sorted(f(1, 2, n)) def a(n): if n == 1: return 1 x = primorial(n) y = 2*x while True: v = omega_cond(x, y, n) if len(v) > 0: return integer_nthroot(v[0], 2)[0] print("done to y =", y, time()-time0) x = y+1 y = 2*x print([a(n) for n in range(1, 5)]) from time import time time0 = time() alst = [] for n in range(1, 200): an = a(n) alst.append(an) print(n, an, len(str(alst))-2, time()-time0) print(" ", alst) print(" ", data)