jordan_polya_numbers.sf

#!/usr/bin/ruby

# Jordan-Pólya numbers (or J-P numbers for short) are the numbers that can be obtained by multiplying together one or more (not necessarily distinct) factorials.

# See also:
#   https://oeis.org/A001013
#   https://rosettacode.org/wiki/Jordan-P%C3%B3lya_numbers
#   https://en.wikipedia.org/wiki/Jordan%E2%80%93P%C3%B3lya_number

func jordan_pólya_numbers(n, k, F, callback) {

    var factors = F.sort.uniq
    var factors_end = factors.end

    if (k == 0) {
        callback(1)
        return nil
    }

    func (m, k, i=0) {

        if (k == 1) {

            var L = idiv(n,m)

            for j in (i..factors_end) {
                with (factors[j]) {|q|
                    q > L && break
                    callback(m*q)
                }
            }

            return nil
        }

        var L = idiv(n,m).iroot(k)

        for j in (i..factors_end) {
            with (factors[j]) { |q|
                q > L && break
                __FUNC__(m*q, k-1, j)
            }
        }
    }(1, k)

    return nil
}

func inverse_factorial_W(n) {
    var l = (log(n) - log(Num.tau)/2)
    l / lambert_w(l / Num.e) - 1/2
}

var limit = 1e9
var factors = (2..inverse_factorial_W(limit).int -> map { .factorial })
var terms = Set()

for k in (0 .. limit.ilog2) {
    jordan_pólya_numbers(limit, k, factors, {|v| terms << v })
}

terms.sort!

say "The first 50 Jordan-Pólya numbers:"
terms.first(50).each_slice(10, {|*a|
    a.map { '%5s' % _ }.join(' ').say
})

say "\nThere are #{terms.len} Jordan-Pólya numbers <= #{limit.commify}, where largest is #{terms.last}."

__END__
The first 50 Jordan-Pólya numbers:
    1     2     4     6     8    12    16    24    32    36
   48    64    72    96   120   128   144   192   216   240
  256   288   384   432   480   512   576   720   768   864
  960  1024  1152  1296  1440  1536  1728  1920  2048  2304
 2592  2880  3072  3456  3840  4096  4320  4608  5040  5184

There are 528 Jordan-Pólya numbers <= 1,000,000,000, where largest is 995328000.