# Python program for OEIS A218030 # Michael S. Branicky, Oct 06 2022 # A218030 Numbers k equal to half of the product of the nonzero (base-10) digits of k^2. +40 data = [2, 5, 54, 648, 2160, 337169025526136832000, 685506275314921762068267522458966662115416623590907309075726336000000, 46641846972427276691124922228108091690332947069125333309512419901440000000000] import heapq from sympy import primerange, integer_nthroot from time import time time0 = time() def nzprod(n): p = 1 for i in map(int, (d for d in str(n) if d != "0")): p *= i return p def psmoothgen(p, m=1): # generate all p-smooth terms multiples of m v, oldv, h, psmooth_primes = 1, 0, [1, p], list(primerange(1, p+1)) while True: v = heapq.heappop(h) if v != oldv: yield v*m oldv = v for p in psmooth_primes: heapq.heappush(h, v*p) def ok(t): return nzprod(t**2) == 2*t def afind(): alst = [] digits = 0 for s in psmoothgen(7, m=2): t = s//2 d = len(str(t)) if d > digits: digits = d print("...", digits, time()-time0) if ok(t): alst.append(t) print("FOUND", t, len(str(alst))-2) print(" ", alst) print(" ", data) with open('b218030.txt', 'a') as bfile: bfile.write(f"{len(alst)} {t}\n") afind()