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.

Jordan-Pólya numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Example

480 is a J-P number because 480 = 2! x 2! x 5!.

Task

Find and show on this page the first 50 J-P numbers.

What is the largest J-P number less than 100 million?

Bonus

Find and show on this page the 800th, 1,800th, 2,800th and 3,800th J-P numbers and also show their decomposition into factorials in highest to lowest order.

Hint: These J-P numbers are all less than 2^53.

References


jq

Works with: jq

Also works with gojq, the Go implementation of jq

Adapted from Wren

### Generic functions
# For gojq
def _nwise($n):
  def n: if length <= $n then . else .[0:$n] , (.[$n:] | n) end;
  n;

def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;

# tabular print
def tprint(columns; wide):
  reduce _nwise(columns) as $row ("";
     . + ($row|map(lpad(wide)) | join(" ")) + "\n" );

# Input: an array
# Output: a stream of pairs [$x, $frequency]
# A two-level dictionary is used: .[type][tostring]
def frequencies:
  if length == 0 then empty
  else . as $in
  | reduce range(0; length) as $i ({};
     $in[$i] as $x
     | .[$x|type][$x|tostring] as $pair
     | if $pair
       then .[$x|type][$x|tostring] |= (.[1] += 1)
       else .[$x|type][$x|tostring] = [$x, 1]
       end )
  | .[][]
  end ;

# Output: the items in the stream up to but excluding the first for which cond is truthy
def emit_until(cond; stream): label $out | stream | if cond then break $out else . end;

### Jordan-Pólya numbers
# input: {factorial}
# output: an array
def JordanPolya($lim; $mx):
  if $lim < 2 then [1]
  else . + {v: [1], t: 1, k: 2}
  | .mx = ($mx // $lim)
  | until(.k > .mx or .t > $lim;
        .t *= .k
	| if .t <= $lim
          then reduce JordanPolya(($lim/.t)|floor; .t)[] as $rest (.;
                 .v += [.t * $rest] ) 
          | .k += 1
	  else .
	  end)
  | .v	
  | unique
  end;

# Cache m! for m <= $n
def cacheFactorials($n):
   {fact: 1, factorial: [1]}
   | reduce range(1; $n + 1) as $i (.;
       .fact *= $i
       | .factorial[$i] = .fact );

# input: {factorial}
def Decompose($n; $start):
  if $start and $start < 2 then []
  else 
  { factorial,
    start: ($start // 18),
    m: $n,
    f: [] }
  | label $out    
  | foreach range(.start; 1; -1) as $i (.;
        .i = $i
        | .emit = null
        | until (.emit or (.m % .factorial[$i] != 0);
            .f += [$i]
            | .m = (.m / .factorial[$i])
            | if .m == 1 then .emit = .f else . end)
	| if .emit then ., break $out else . end)
  | if .emit then .emit
    elif .i == 2 then Decompose($n; .start-1)
    else empty
    end
  end;

# Input: {factorial}
# $v should be an array of J-P numbers
def synopsis($v):
  (100, 800, 1800, 2800, 3800) as $i
  | if $v[$i-1] == null 
    then "\nThe \($i)th Jordan-Pólya number was not found." | error
    else "\nThe \($i)th Jordan-Pólya number is \($v[$i-1] )",
          ([Decompose($v[$i-1]; null) | frequencies]
           | map( if (.[1] == 1) then "\(.[0])!"  else "(\(.[0])!)^\(.[1])" end)
           | "  i.e. " + join(" * ") )
    end ;

def task:
  cacheFactorials(18)
  | JordanPolya(pow(2;53)-1; null) as $v
  | "\($v|length) Jordan–Pólya numbers have been found. The first 50 are:",
    ( $v[:50] | tprint(10; 4)),
    "\nThe largest Jordan–Pólya number before 100 million: " +
    "\(if $v[-1] > 1e8 then last(emit_until(. >= 1e8; $v[])) else "not found" end)",
    synopsis($v) ;

task
Output:

gojq and jq produce the same results except that gojq produces the factorizations in a different order. The output shown here corresponds to the invocation: jq -nr -f jordan-polya-numbers.jq

3887 Jordan–Pólya numbers have been found. The first 50 are:
   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


The largest Jordan–Pólya number before 100 million: 99532800

The 100th Jordan-Pólya number is 92160
  i.e. 6! * (2!)^7

The 800th Jordan-Pólya number is 18345885696
  i.e. (4!)^7 * (2!)^2

The 1800th Jordan-Pólya number is 9784472371200
  i.e. (6!)^2 * (4!)^2 * (2!)^15

The 2800th Jordan-Pólya number is 439378587648000
  i.e. 14! * 7!

The 3800th Jordan-Pólya number is 7213895789838336
  i.e. (4!)^8 * (2!)^16

Julia

""" The aupto() function is taken from the Python code at oeis.org/A001013 """

function aupto(lim::T, mx::T = zero(T)) where T <: Integer
    lim < 2 && return [one(T)]
    v, t = [one(T)], one(T)
    mx == 0 && (mx = lim)
    for k in 2:mx
        t *= k
        t > lim && break
        append!(v, [t * rest for rest in aupto(lim ÷ t, t)])
    end
    return unique(sort!(v))
end

const factorials = map(factorial, 2:18)

function factor_as_factorials(n::T) where T <: Integer
    fac_exp = Tuple{Int, Int}[]
    for idx in length(factorials):-1:1
        m = n
        empty!(fac_exp)
        for i in idx:-1:1
            expo = 0
            while m % factorials[i] == 0
                expo += 1
                m ÷= factorials[i]
            end
            if expo > 0
                push!(fac_exp, (i + 1, expo))
            end
        end
        m == 1 && break
    end
    return fac_exp
end

const superchars = ["\u2070", "\u00b9", "\u00b2", "\u00b3", "\u2074",
                    "\u2075", "\u2076", "\u2077", "\u2078", "\u2079"]
super(n::Integer) = prod([superchars[i + 1] for i in reverse(digits(n))])

arr = aupto(2^53)

println("First 50 Jordan–Pólya numbers:")
foreach(p -> print(rpad(p[2], 6), p[1] % 10 == 0 ? "\n" : ""), enumerate(arr[1:50]))

println("\nThe largest Jordan–Pólya number before 100 million: ", arr[findlast(<(100_000_000), arr)])

for n in [800, 1800, 2800, 3800]
    print("\nThe $(n)th Jordan-Pólya number is: $(arr[n])\n= ")
    println(join(map(t -> "($(t[1])!)$(super(t[2])))", factor_as_factorials(arr[n])), " x "))
end
Output:
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

The largest Jordan–Pólya number before 100 million: 99532800

The 800th Jordan-Pólya number is: 18345885696
= (4!)⁷) x (2!)²)

The 1800th Jordan-Pólya number is: 9784472371200
= (6!)²) x (4!)²) x (2!)¹⁵)

The 2800th Jordan-Pólya number is: 439378587648000
= (14!)¹) x (7!)¹)

The 3800th Jordan-Pólya number is: 7213895789838336
= (4!)⁸) x (2!)¹⁶)

Nim

import std/[algorithm, math, sequtils, strformat, strutils, tables]

const Max = if sizeof(int) == 8: 20 else: 12

type Decomposition = CountTable[int]

const Superscripts: array['0'..'9', string] = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]

func superscript(n: Natural): string =
  ## Return the Unicode string to use to represent an exponent.
  if n == 1:
    return ""
  for d in $n:
    result.add Superscripts[d]

proc `$`(d: Decomposition): string =
  ## Return the representation of a decomposition.
  for (value, count) in sorted(d.pairs.toSeq, Descending):
    result.add &"({value}!){superscript(count)}"


# List of Jordan-Pólya numbers and their decomposition.
var jordanPolya = @[1]
var decomposition: Table[int, CountTable[int]] = {1: initCountTable[int]()}.toTable

# Build the list and the decompositions.
for m in 2..Max:                  # Loop on each factorial.
  let f = fac(m)
  for k in 0..jordanPolya.high:   # Loop on existing elements.
    var n = jordanPolya[k]
    while n <= int.high div f:    # Multiply by successive powers of n!
      let p = n
      n *= f
      jordanPolya.add n
      decomposition[n] = decomposition[p]
      decomposition[n].inc(m)

# Sort the numbers and remove duplicates.
jordanPolya = sorted(jordanPolya).deduplicate(true)

echo "First 50 Jordan-Pólya numbers:"
for i in 0..<50:
  stdout.write &"{jordanPolya[i]:>4}"
  stdout.write if i mod 10 == 9: '\n' else: ' '

echo "\nLargest Jordan-Pólya number less than 100 million: ",
     insertSep($jordanPolya[jordanPolya.upperBound(100_000_000) - 1])

for i in [800, 1800, 2800, 3800]:
  let n = jordanPolya[i - 1]
  var d = decomposition[n]
  echo &"\nThe {i}th Jordan-Pólya number is:"
  echo &"{insertSep($n)} = {d}"
Output:
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

Largest Jordan-Pólya number less than 100 million: 99_532_800

The 800th Jordan-Pólya number is:
18_345_885_696 = (4!)⁷(2!)²

The 1800th Jordan-Pólya number is:
9_784_472_371_200 = (6!)²(4!)²(2!)¹⁵

The 2800th Jordan-Pólya number is:
439_378_587_648_000 = (14!)(7!)

The 3800th Jordan-Pólya number is:
7_213_895_789_838_336 = (4!)⁸(2!)¹⁶

Phix

with javascript_semantics
function factorials_le(atom limit)
    sequence res = {}
    while true do
        atom nf = factorial(length(res)+1)
        if nf>limit then exit end if
        res &= nf
    end while
    return res
end function

function jp(atom limit)
    sequence res = factorials_le(limit)
    integer k=2
    while k<=length(res) do
        atom rk = res[k]
        for l=2 to length(res) do
            atom kl = res[l]*rk
            if kl>limit then exit end if
            do
                integer p = binary_search(kl,res)
                if p<0 then
                    p = abs(p)
                    res = res[1..p-1] & kl & res[p..$]
                end if
                kl *= rk
            until kl>limit
        end for
        k += 1
    end while
    return res
end function

function decompose(atom jp)
    atom jp0 = jp
    sequence f = factorials_le(jp)
    while true do
        string res = ""
        for i=length(f) to 2 by -1 do
            atom fi = f[i]
            if remainder(jp,fi)=0 then
                jp /= fi
                if length(res) then res &= " * " end if
                res &= sprintf("%d!",i)
                integer fc = 0
                while remainder(jp,fi)=0 do
                    jp /= fi
                    fc += 1
                end while
                if fc then
                    res &= sprintf("^%d",fc+1)
                end if
            end if
            if jp=1 then return res end if
        end for
        f = f[1..-2]
        jp = jp0
    end while
end function

atom t0 = time()
sequence r = jp(power(2,53)-1)
printf(1,"%d Jordan-Polya numbers found, the first 50 are:\n%s\n",
         {length(r),join_by(r[1..50],1,10," ",fmt:="%4d")})
printf(1,"The largest under 100 million: %,d\n",r[abs(binary_search(1e8,r))-1])
for i in {100,800,1800,2800,3800} do
    printf(1,"The %d%s is %,d = %s\n",{i,ord(i),r[i],decompose(r[i])})
end for
?elapsed(time()-t0)
Output:
3887 Jordan-Polya numbers found, the first 50 are:
   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

The largest under 100 million: 99,532,800
The 100th is 92,160 = 6! * 2!^7
The 800th is 18,345,885,696 = 4!^7 * 2!^2
The 1800th is 9,784,472,371,200 = 6!^2 * 4!^2 * 2!^15
The 2800th is 439,378,587,648,000 = 14! * 7!
The 3800th is 7,213,895,789,838,336 = 4!^8 * 2!^16
"0.3s"

Actually slightly faster under pwa/p2js than it is on desktop/Phix.

Wren

Library: Wren-set
Library: Wren-seq
Library: Wren-fmt

This uses the recursive PARI/Python algorithm in the OEIS entry.

import "./set" for Set
import "./seq" for Lst
import "./fmt" for Fmt

var JordanPolya = Fn.new { |lim, mx|
    if (lim < 2) return [1]
    var v = Set.new()
    v.add(1)
    var t = 1
    if (!mx) mx = lim
    for (k in 2..mx) {
        t = t * k
        if (t > lim) break
        for (rest in JordanPolya.call((lim/t).floor, t)) {
            v.add(t * rest)
        }
    }
    return v.toList.sort()
}

var factorials = List.filled(19, 1)

var cacheFactorials = Fn.new {
    var fact = 1
    for (i in 2..18) {
        fact = fact * i
        factorials[i] = fact
    }
}

var Decompose = Fn.new { |n, start|
    if (!start) start = 18
    if (start < 2) return []
    var m = n
    var f = []
    for (i in start..2) {
        while (m % factorials[i] == 0) {
            f.add(i)
            m =  m / factorials[i]
            if (m == 1) return f
        }
    }
    return Decompose.call(n, start-1)
}

cacheFactorials.call()
var v = JordanPolya.call(2.pow(53)-1, null)
System.print("First 50 Jordan–Pólya numbers:")
Fmt.tprint("$4d ", v[0..49], 10)

System.write("\nThe largest Jordan–Pólya number before 100 million: ")
for (i in 1...v.count) {
    if (v[i] > 1e8) {
        Fmt.print("$,d\n", v[i-1])
        break
    }
}

for (i in [800, 1800, 2800, 3800]) {
    Fmt.print("The $,r Jordan-Pólya number is : $,d", i, v[i-1])
    var g = Lst.individuals(Decompose.call(v[i-1], null))
    var s = g.map { |l|
        if (l[1] == 1) return "%(l[0])!"
        return Fmt.swrite("($d!)$S", l[0], l[1])
    }.join(" x ")
    Fmt.print("= $s\n", s)
}
Output:
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 

The largest Jordan–Pólya number before 100 million: 99,532,800

The 800th Jordan-Pólya number is : 18,345,885,696
= (4!)⁷ x (2!)²

The 1,800th Jordan-Pólya number is : 9,784,472,371,200
= (6!)² x (4!)² x (2!)¹⁵

The 2,800th Jordan-Pólya number is : 439,378,587,648,000
= 14! x 7!

The 3,800th Jordan-Pólya number is : 7,213,895,789,838,336
= (4!)⁸ x (2!)¹⁶

XPL0

Simple-minded brute force. 20 seconds on Pi4. No bonus.

int Factorials(1+12);

func IsJPNum(N0);
int  N0;
int  N, Limit, I, Q;
[Limit:= 7;
N:= N0;
loop    [I:= Limit;
        loop    [Q:= N / Factorials(I);
                if rem(0) = 0 then
                        [if Q = 1 then return true;
                        N:= Q;
                        ]
                else    I:= I-1;
                if I = 1 then 
                        [if Limit = 1 then return false;
                        Limit:= Limit-1;
                        N:= N0;
                        quit;
                        ]
                ];
        ];
];

int F, N, C, SN;
[F:= 1;
for N:= 1 to 12 do
        [F:= F*N;
        Factorials(N):= F;
        ];
Text(0, "First 50 Jordan-Polya numbers:^m^j");
Format(5, 0);
RlOut(0, 1.);           \handle odd number exception
C:= 1;  N:= 2;  
loop    [if IsJPNum(N) then
                [C:= C+1;
                if C <= 50 then
                        [RlOut(0, float(N));
                        if rem(C/10) = 0 then CrLf(0);
                        ];
                SN:= N;
                ];
        N:= N+2;
        if N >= 100_000_000 then quit;
        ];
Text(0, "^m^jThe largest Jordan-Polya number before 100 million: ");
IntOut(0, SN);  CrLf(0);
]
Output:
First 50 Jordan-Polya 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

The largest Jordan-Polya number before 100 million: 99532800