# 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()