(PARI) is(n) = { if(n == 1, return(0)); my(i, factorials, e, res, v); f = factor(n); if(prime(#f~) != f[#f~, 1], return(0); ); if(f[, 2] != vecsort(f[, 2], , 4), return(0); ); factorials = List(); e = List(); res = List(); for(i = 2, oo, v = valuation(n, i!); if(v > 0, listput(factorials, i!); listput(e, v); , break ) ); forvec(x = vector(#e-1, i, [0, e[i+1]]), c = prod(i = 1, #e-1, factorials[i+1]^x[i]); if(c <= n && denominator(n/c) == 1&& 1 << logint(n/c, 2) == n/c, listput(res, concat([valuation(n/c, 2)], x)) ) ); #res >= 2 } \\
David A. Corneth, Jan 13 2023
(Python)
from itertools import count, islice
from functools import lru_cache
from sympy import integer_log
from oeis_sequences.OEISsequences import bisection
def
A359750_gen(): # generator of terms
@lru_cache(maxsize=None)
def g(x, m, j): return sum(g(x//m**i, m-1, i) for i in range(j, integer_log(x, m)[0]+1)) if m-2 else max(0, x.bit_length()-j)
def f(x):
c, f = n-1+x, 1
for k in count(2):
f *= k
if f>x:
break
c -= g(x, k, 1)
return c
a, b = 0, 0
for n in count(1):
k = bisection(f, max(b, n), max(b, n))
if k == b and k != a:
yield k
a, b = b, k
