Jordan-Pólya numbers

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.

Task
Jordan-Pólya numbers
You are encouraged to solve this task according to the task description, using any language you may know.
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. Optionally, do the same for the 1,050th J-P number.

Where there is more than one way to decompose a J-P number into factorials, choose the way which uses the largest factorials.

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

References


C

Translation of: Wren
Library: GLib

A translation of the second version. Run time about 0.035 seconds.

#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <locale.h>
#include <glib.h>

uint64_t factorials[19] = {1, 1};

void cacheFactorials() {
    uint64_t fact = 1;
    int i;
    for (i = 2; i < 19; ++i) {
        fact *= i;
        factorials[i] = fact;
    }
}

int findNearestFact(uint64_t n) {
    int i;
    for (i = 1; i < 19; ++i) {
        if (factorials[i] >= n) return i;
    }
    return 18;
}

int findNearestInArray(GArray *a, uint64_t n) {
    int l = 0, r = a->len, m;
    while (l < r) {
        m = (l + r)/2;
        if (g_array_index(a, uint64_t, m) > n) {
            r = m;
        } else {
            l = m + 1;
        }
    }
    if (r > 0 && g_array_index(a, uint64_t, r-1) == n) return r - 1;
    return r;
}

GArray *jordanPolya(uint64_t limit) {
    int i, ix, k, l, p;
    uint64_t t, rk, kl;
    GArray *res = g_array_new(false, false, sizeof(uint64_t));
    ix = findNearestFact(limit);
    for (i = 0; i <= ix; ++i) {
        t = factorials[i];
        g_array_append_val(res, t);
    }
    k = 2;
    while (k < res->len) {
        rk = g_array_index(res, uint64_t, k);
        for (l = 2; l < res->len; ++l) {
            t = g_array_index(res, uint64_t, l);
            if (t > limit/rk) break;
            kl = t * rk;
            while (true) {
                p = findNearestInArray(res, kl);
                if (p < res->len && g_array_index(res, uint64_t, p) != kl) {
                    g_array_insert_val(res, p, kl);
                } else if (p == res->len) {
                    g_array_append_val(res, kl);
                }
                if (kl > limit/rk) break;
                kl *= rk;
            }
        }
        ++k;
    }
    return g_array_remove_index(res, 0);
}

GArray *decompose(uint64_t n, int start) {
    uint64_t i, s, t, m, prod;
    GArray *f, *g;
    for (s = start; s > 0; --s) {
        f = g_array_new(false, false, sizeof(uint64_t));
        if (s < 2) return f;
        m = n;
        while (!(m % factorials[s])) {
            g_array_append_val(f, s);
            m /= factorials[s];
            if (m == 1) return f;
        }
        if (f->len > 0) {
            g = decompose(m, s - 1);
            if (g->len > 0) {
                prod = 1;
                for (i = 0; i < g->len; ++i) {
                    prod *= factorials[(int)g_array_index(g, uint64_t, i)];
                }
                if (prod == m) {
                    for (i = 0; i < g->len; ++i) {
                        t = g_array_index(g, uint64_t, i);
                        g_array_append_val(f, t);
                    }
                    g_array_free(g, true);
                    return f;
                }
            }
            g_array_free(g, true);
        }
        g_array_free(f, true);
    }
}

char *superscript(int n) {
    char* ss[] = {"⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"};
    if (n < 10) return ss[n];
    static char buf[7];
    sprintf(buf, "%s%s", ss[n/10], ss[n%10]);
    return buf;
}
                 
int main() {
    int i, j, ix, count;
    uint64_t t, u;
    GArray *v, *w;
    cacheFactorials();
    v = jordanPolya(1ull << 53);
    printf("First 50 Jordan-Pólya numbers:\n");
    for (i = 0; i < 50; ++i) {
        printf("%4ju ", g_array_index(v, uint64_t, i));
        if (!((i + 1) % 10)) printf("\n");
    }
    printf("\nThe largest Jordan-Pólya number before 100 millon: ");
    setlocale(LC_NUMERIC, "");
    ix = findNearestInArray(v, 100000000ull);    
    printf("%'ju\n\n", g_array_index(v, uint64_t, ix - 1));

    uint64_t targets[5] = {800, 1050, 1800, 2800, 3800};
    for (i = 0; i < 5; ++i) {
        t = g_array_index(v, uint64_t, targets[i] - 1);
        printf("The %'juth Jordan-Pólya number is : %'ju\n", targets[i], t);
        w = decompose(t, 18);
        count = 1;
        t = g_array_index(w, uint64_t, 0);
        printf(" = ");
        for (j = 1; j < w->len; ++j) {
            u = g_array_index(w, uint64_t, j);
            if (u != t) {
                if (count == 1) {
                    printf("%ju! x ", t);
                } else {
                    printf("(%ju!)%s x ", t, superscript(count));
                    count = 1;
                }
                t = u;
            } else {
                ++count;
            }
        }
        if (count == 1) {
            printf("%ju! x ", t);
        } else {
            printf("(%ju!)%s x ", t, superscript(count));
        }
        printf("\b\b \n\n");
        g_array_free(w, true);
    }
    g_array_free(v, true);
    return 0;
}
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 millon: 99,532,800

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

The 1,050th Jordan-Pólya number is : 139,345,920,000
 = 8! x (5!)³ 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!)¹⁶   

C#

Translation of: Java
using System;
using System.Collections.Generic;
using System.Linq;

public class JordanPolyaNumbers
{
    private static SortedSet<long> jordanPolyaSet = new SortedSet<long>();
    private static Dictionary<long, SortedDictionary<int, int>> decompositions = new Dictionary<long, SortedDictionary<int, int>>();

    public static void Main(string[] args)
    {
        CreateJordanPolya();

        long belowHundredMillion = jordanPolyaSet.LastOrDefault(x => x < 100_000_000L);
        List<long> jordanPolya = new List<long>(jordanPolyaSet);

        Console.WriteLine("The first 50 Jordan-Polya numbers:");
        for (int i = 0; i < 50; i++)
        {
            Console.Write($"{jordanPolya[i],5}{(i % 10 == 9 ? "\n" : "")}");
        }
        Console.WriteLine();

        Console.WriteLine("The largest Jordan-Polya number less than 100 million: " + belowHundredMillion);
        Console.WriteLine();

        foreach (int i in new List<int> { 800, 1050, 1800, 2800, 3800 })
        {
            Console.WriteLine($"The {i}th Jordan-Polya number is: {jordanPolya[i - 1]} = {ToString(decompositions[jordanPolya[i - 1]])}");
        }
    }

    private static void CreateJordanPolya()
    {
        jordanPolyaSet.Add(1L);
        SortedSet<long> nextSet = new SortedSet<long>();
        decompositions[1L] = new SortedDictionary<int, int>();
        long factorial = 1;

        for (int multiplier = 2; multiplier <= 20; multiplier++)
        {
            factorial *= multiplier;
            foreach (long number in new SortedSet<long>(jordanPolyaSet))
            {
                long newNumber = number;
                while (newNumber <= long.MaxValue / factorial)
                {
                    long original = newNumber;
                    newNumber *= factorial;
                    nextSet.Add(newNumber);

                    decompositions[newNumber] = new SortedDictionary<int, int>(decompositions[original]);
                    if (decompositions[newNumber].ContainsKey(multiplier))
                    {
                        decompositions[newNumber][multiplier]++;
                    }
                    else
                    {
                        decompositions[newNumber][multiplier] = 1;
                    }
                }
            }
            jordanPolyaSet.UnionWith(nextSet);
            nextSet.Clear();
        }
    }

    private static string ToString(SortedDictionary<int, int> map)
    {
        string result = "";
        foreach (int key in map.Keys)
        {
            result = key + "!" + (map[key] == 1 ? "" : "^" + map[key]) + " * " + result;
        }
        return result.TrimEnd(' ', '*');
    }
}
Output:
The 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 less than 100 million: 99532800

The 800th Jordan-Polya number is: 18345885696 = 4!^7 * 2!^2
The 1050th Jordan-Polya number is: 139345920000 = 8! * 5!^3 * 2!
The 1800th Jordan-Polya number is: 9784472371200 = 6!^2 * 4!^2 * 2!^15
The 2800th Jordan-Polya number is: 439378587648000 = 14! * 7!
The 3800th Jordan-Polya number is: 7213895789838336 = 4!^8 * 2!^16

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <map>
#include <set>
#include <unordered_map>
#include <vector>

constexpr int64_t LIMIT = static_cast<uint64_t>(1) << 53;

std::set<int64_t> jordan_polya_set;
std::unordered_map<int64_t, std::map<int32_t, int32_t>> decompositions;

std::string toString(const std::map<int32_t, int32_t>& a_map) {
	std::string result;
	for ( const auto& [key, value] : a_map ) {
		result = std::to_string(key) + "!" + ( value == 1 ? "" : "^" + std::to_string(value) ) + " * " + result;
	}
	return result.substr(0, result.length() - 3);
}

std::vector<int64_t> set_to_vector(const std::set<int64_t>& a_set) {
    std::vector<int64_t> result;
    result.reserve(a_set.size());

    for ( const int64_t& element : a_set ) {
        result.emplace_back(element);
    }
    return result;
}

void insert_or_update(std::map<int32_t, int32_t>& map, const int32_t& entry) {
	if ( map.find(entry) == map.end() ) {
		map.emplace(entry, 1);
	} else {
		map[entry]++;
	}
}

void create_jordan_polya() {
	jordan_polya_set.emplace(1);
	decompositions[1] = std::map<int32_t, int32_t>();
	int64_t factorial = 1;

	for ( int32_t multiplier = 2; multiplier <= 20; ++multiplier ) {
		factorial *= multiplier;
		for ( int64_t number : jordan_polya_set ) {
			while ( number <= LIMIT / factorial ) {
				int64_t original = number;
				number *= factorial;
				jordan_polya_set.emplace(number);
				decompositions[number] = decompositions[original];
				insert_or_update(decompositions[number], multiplier);
			}
		}
	}
}

int main() {
	create_jordan_polya();

	std::vector<int64_t> jordan_polya = set_to_vector(jordan_polya_set);

	std::cout << "The first 50 Jordan-Polya numbers:" << std::endl;
	for ( int64_t i = 0; i < 50; ++i ) {
		std::cout << std::setw(5) << jordan_polya[i] << ( i % 10 == 9 ? "\n" : "" );
	}

	const std::vector<int64_t>::iterator hundred_million =
		std::lower_bound(jordan_polya.begin(), jordan_polya.end(), 100'000'000);
	std::cout << "\n" << "The largest Jordan-Polya number less than 100 million: "
			  << jordan_polya[(hundred_million - jordan_polya.begin() - 1)] << std::endl << std::endl;

	for ( int32_t i : { 800, 1050, 1800, 2800, 3800 } ) {
		std::cout << "The " << i << "th Jordan-Polya number is: " << jordan_polya[i - 1]
			      << " = " << toString(decompositions[jordan_polya[i - 1]]) << std::endl;
	}
}
Output:
The 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 less than 100 million: 99532800

The 800th Jordan-Polya number is: 18345885696 = 4!^7 * 2!^2
The 1050th Jordan-Polya number is: 139345920000 = 8! * 5!^3 * 2!
The 1800th Jordan-Polya number is: 9784472371200 = 6!^2 * 4!^2 * 2!^15
The 2800th Jordan-Polya number is: 439378587648000 = 14! * 7!
The 3800th Jordan-Polya number is: 7213895789838336 = 4!^8 * 2!^16

Dart

Translation of: Java
import "dart:collection";
import "dart:io";

void main() {
  createJordanPolya();

  final belowHundredMillion = jordanPolyaSet.lastWhere((element) => element <= 100000000, orElse: () => null);
  List<int> jordanPolya = jordanPolyaSet.toList();

  print("The first 50 Jordan-Polya numbers:");
  for (int i = 0; i < 50; i++) {
    // Right-align each number in a 5-character wide space
    stdout.write(jordanPolya[i].toString().padLeft(5));
    if (i % 10 == 9) {
      print(""); // Newline every 10 numbers
    }
  }
  print("");

  print("The largest Jordan-Polya number less than 100 million: ${belowHundredMillion ?? 'Not found'}");
  print("");

  for (int i in [800, 1050, 1800, 2800, 3800]) {
    var decomposition = decompositions[jordanPolya[i - 1]];
    if (decomposition != null) {
      print("The ${i}th Jordan-Polya number is: ${jordanPolya[i - 1]} = ${mapToString(decomposition)}");
    }
  }
}

SplayTreeSet<int> jordanPolyaSet = SplayTreeSet<int>();
Map<int, Map<int, int>> decompositions = {};

void createJordanPolya() {
  jordanPolyaSet.add(1);
  Set<int> nextSet = SplayTreeSet<int>();
  decompositions[1] = {};
  int factorial = 1;

  for (int multiplier = 2; multiplier <= 20; multiplier++) {
    factorial *= multiplier;
    for (int number in jordanPolyaSet) {
      int tempNumber = number;
      while (tempNumber <= 9223372036854775807 ~/ factorial) {
        int original = tempNumber;
        tempNumber *= factorial;
        nextSet.add(tempNumber);
        var originalDecomposition = decompositions[original];
        if (originalDecomposition != null) {
          decompositions[tempNumber] = Map<int, int>.from(originalDecomposition)
            ..update(multiplier, (value) => value + 1, ifAbsent: () => 1);
        }
      }
    }
    jordanPolyaSet.addAll(nextSet);
    nextSet.clear();
  }
}

String mapToString(Map<int, int> map) {
  String result = "";
  map.forEach((key, value) {
    result = "$key!${value == 1 ? '' : '^$value'} * $result";
  });
  return result.isEmpty ? result : result.substring(0, result.length - 3);
}
Output:
The 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 less than 100 million: 99532800

The 800th Jordan-Polya number is: 18345885696 = 4!^7 * 2!^2
The 1050th Jordan-Polya number is: 139345920000 = 8! * 5!^3 * 2!
The 1800th Jordan-Polya number is: 9784472371200 = 6!^2 * 4!^2 * 2!^15
The 2800th Jordan-Polya number is: 439378587648000 = 14! * 7!
The 3800th Jordan-Polya number is: 7213895789838336 = 4!^8 * 2!^16

EasyLang

Translation of: FreeBASIC
fastfunc jpnum m .
   n = m
   limite = 7
   while 1 = 1
      fac = 1
      i = 1
      while i < limite
         i += 1
         fac *= i
      .
      repeat
         q = n div fac
         if n mod fac = 0
            if q = 1
               return 1
            .
            n = q
         else
            fac = fac / i
            i -= 1
         .
         until i = 1
      .
      limite -= 1
      if limite = 0
         return 0
      .
      n = m
   .
.
numfmt 0 5
write 1
c = 1
n = 2
repeat
   if jpnum n = 1
      c += 1
      if c <= 50
         write n
         if c mod 8 = 0
            print ""
         .
      .
      sn = n
   .
   n += 2
   until n >= 1e8
.
print ""
print "The largest Jordan-Polya number before 100 million: " & sn

FreeBASIC

Translation of: XPL0

Simple-minded brute force. No bonus.

Dim Shared As Uinteger Factorials(1+12)

Function isJPNum(m As Uinteger) As Boolean
    Dim As Uinteger n = m, limite = 7, i, q
    Do
        i = limite
        Do
            q = n / Factorials(i)
            If n Mod Factorials(i) = 0 Then
                If q = 1 Then Return True
                n = q
            Else    
                i -= 1
            End If
            If i = 1 Then 
                If limite = 1 Then Return False
                limite -= 1
                n = m
                Exit Do
            End If
        Loop
    Loop
End Function

Dim As Uinteger fact = 1, n
For n = 1 To 12
    fact *= n
    Factorials(n) = fact
Next

Print "First 50 Jordan-Polya numbers:"
Print "    1";
Dim As Uinteger c, sn
c = 1  
n = 2 
Do
    If isJPNum(n) Then
        c += 1
        If c <= 50 Then
            Print Using "#####"; n;
            If c Mod 10 = 0 Then Print
        End If
        sn = n
    End If
    n += 2
Loop Until n >= 1e8

Print !"\nThe largest Jordan-Polya number before 100 million: "; sn

Sleep
Output:
Same as XPL0 entry.

Go

Translation of: C
Library: Go-rcu

Run time a shade slower than C at about 0.047 seconds.

package main

import (
    "fmt"
    "rcu"
    "sort"
)

var factorials = make([]uint64, 19)

func cacheFactorials() {
    factorials[0] = 1
    for i := uint64(1); i < 19; i++ {
        factorials[i] = factorials[i-1] * i
    }
}

func jordan_polya(limit uint64) []uint64 {
    ix := sort.Search(19, func(i int) bool { return factorials[i] >= limit })
    if ix > 18 {
        ix = 18
    }
    var res []uint64
    res = append(res, factorials[0:ix+1]...)
    k := 2
    for k < len(res) {
        rk := res[k]
        for l := 2; l < len(res); l++ {
            t := res[l]
            if t > limit/rk {
                break
            }
            kl := t * rk
            for {
                p := sort.Search(len(res), func(i int) bool { return res[i] >= kl })
                if p < len(res) && res[p] != kl {
                    res = append(res[0:p+1], res[p:]...)
                    res[p] = kl
                } else if p == len(res) {
                    res = append(res, kl)
                }
                if kl > limit/rk {
                    break
                }
                kl *= rk
            }
        }
        k++
    }
    return res[1:]
}

func decompose(n uint64, start int) []uint64 {
    for s := uint64(start); s > 0; s-- {
        var f []uint64
        if s < 2 {
            return f
        }
        m := n
        for m%factorials[s] == 0 {
            f = append(f, s)
            m /= factorials[s]
            if m == 1 {
                return f
            }
        }
        if len(f) > 0 {
            g := decompose(m, int(s-1))
            if len(g) > 0 {
                prod := uint64(1)
                for _, e := range g {
                    prod *= factorials[e]
                }
                if prod == m {
                    return append(f, g...)
                }
            }
        }
    }
    return []uint64{}
}

func superscript(n int) string {
    ss := []string{"⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"}
    if n < 10 {
        return ss[n]
    }
    return ss[n/10] + ss[n%10]
}

func main() {
    cacheFactorials()
    v := jordan_polya(uint64(1) << 53)
    fmt.Println("First 50 Jordan-Pólya numbers:")
    for i := 0; i < 50; i++ {
        fmt.Printf("%4d ", v[i])
        if (i+1)%10 == 0 {
            fmt.Println()
        }
    }
    fmt.Printf("\nThe largest Jordan-Pólya number before 100 million: ")
    ix := sort.Search(len(v), func(i int) bool { return v[i] >= 100_000_000 })
    fmt.Println(rcu.Commatize(v[ix-1]))
    fmt.Println()
    for _, e := range []uint64{800, 1050, 1800, 2800, 3800} {
        fmt.Printf("The %sth Jordan-Pólya number is : %s\n", rcu.Commatize(e), rcu.Commatize(v[e-1]))
        w := decompose(v[e-1], 18)
        count := 1
        t := w[0]
        fmt.Printf(" = ")
        for j := 1; j < len(w); j++ {
            u := w[j]
            if u != t {
                if count == 1 {
                    fmt.Printf("%d! x ", t)
                } else {
                    fmt.Printf("(%d!)%s x ", t, superscript(count))
                    count = 1
                }
                t = u
            } else {
                count++
            }
        }
        if count == 1 {
            fmt.Printf("%d! x ", t)
        } else {
            fmt.Printf("(%d!)%s x ", t, superscript(count))
        }
        fmt.Printf("\b\b \n\n")
    }
}
Output:
Same as C example.

J

F=. !P=. p:i.100x
jpprm=: P{.~F I. 1+]

Fs=. 2}.!i.1+{:P
jpfct=: Fs |.@:{.~ Fs I. 1+]

isjp=: {{
  if. 2>y do. y return.
  elseif. 0 < #(q:y)-.jpprm y do. 0 return.
  else.
    for_f. (#~ ] = <.) (%jpfct) y do.
      if. isjp f do. 1 return. end.
    end.
  end.
  0
}}"0

showjp=: {{
  if. 2>y do. i.0 return. end.
  F=. f{~1 i.~b #inv isjp Y#~b=. (]=<.) Y=. y%f=. jpfct y
  F,showjp y%F
}}

NB. generate a Jordan-Pólya of the given length
jpseq=: {{
  r=. 1 2x   NB. sequence, so far
  f=. 2 6x   NB. factorial factors
  i=. 1 0    NB. index of next item of f for each element of r
  g=. 6 4x   NB. product of r with selected item of f
  while. y>#r do.
    r=. r, nxt=. <./g  NB. next item in r
    j=. I.b=. g=nxt    NB. items of g which just be recalculated 
    if. nxt={:f do.    NB. need new factorial factor/
      f=. f,!2+#f
    end.
    i=. 0,~i+b         NB. update indices into f
    g=. (2*nxt),~((j{r)*((<:#f)<.j{i){f) j} g
  end.
  y{.r
}}

Task:

   5 10$jpseq 50
   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
   <:^:(0=isjp)^:_]1e8
99532800
   showjp 99532800
720 720 24 2 2 2

Note that jp factorizations are not necessarily unique. For example, 120 120 6 6 6 2 2 2 2 2 would also be a jp factorization of 99532800.

Bonus (indicated numbers from jp sequence, followed by a jp factorization):

   s=: jpseq 4000
   (,showjp) (<:800){s
18345885696 24 24 24 24 24 24 24 2 2
   (,showjp) (<:1800){s
9784472371200 720 720 24 24 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
   (,showjp) (<:2800){s
439378587648000 87178291200 5040
   (,showjp) (<:3800){s
7213895789838336 24 24 24 24 24 24 24 24 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.TreeMap;
import java.util.TreeSet;

public final class JordanPolyaNumbers {

	public static void main(String[] aArgs) {		
		createJordanPolya();
		
		final long belowHundredMillion = jordanPolyaSet.floor(100_000_000L);		
		List<Long> jordanPolya = new ArrayList<Long>(jordanPolyaSet);
		
		System.out.println("The first 50 Jordan-Polya numbers:");
		for ( int i = 0; i < 50; i++ ) {
			System.out.print(String.format("%5s%s", jordanPolya.get(i), ( i % 10 == 9 ? "\n" : "" )));
		}
		System.out.println();
		
		System.out.println("The largest Jordan-Polya number less than 100 million: " + belowHundredMillion);
		System.out.println();
		
		for ( int i : List.of( 800, 1050, 1800, 2800, 3800 ) ) {
			System.out.println("The " + i + "th Jordan-Polya number is: " + jordanPolya.get(i - 1)
				+ " = " + toString(decompositions.get(jordanPolya.get(i - 1))));
		}
	}
	
	private static void createJordanPolya() {
		jordanPolyaSet.add(1L);
		Set<Long> nextSet = new TreeSet<Long>();
		decompositions.put(1L, new TreeMap<Integer, Integer>());
		long factorial = 1;
		
		for ( int multiplier = 2; multiplier <= 20; multiplier++ ) {
			factorial *= multiplier;			
			for ( Iterator<Long> iterator = jordanPolyaSet.iterator(); iterator.hasNext(); ) {
				long number = iterator.next();
				while ( number <= Long.MAX_VALUE / factorial ) {
					long original = number;					
					number *= factorial;					
					nextSet.add(number);
		      		decompositions.put(number, new TreeMap<Integer, Integer>(decompositions.get(original)));		      		
		      		decompositions.get(number).merge(multiplier, 1, Integer::sum);
				}
		    }			
			jordanPolyaSet.addAll(nextSet);
			nextSet.clear();
		}
	}
	
	private static String toString(Map<Integer, Integer> aMap) {
		String result = "";
		for ( int key : aMap.keySet() ) {
			result = key + "!" + ( aMap.get(key) == 1 ? "" :"^" + aMap.get(key) ) + " * " + result;
		}		
		return result.substring(0, result.length() - 3);
	}

	private static TreeSet<Long> jordanPolyaSet = new TreeSet<Long>();
	private static Map<Long, Map<Integer, Integer>> decompositions = new HashMap<Long, Map<Integer, Integer>>();	

}
Output:
The 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 less than 100 million: 99532800

The 800th Jordan-Polya number is: 18345885696 = 4!^7 * 2!^2
The 1050th Jordan-Polya number is: 139345920000 = 8! * 5!^3 * 2!
The 1800th Jordan-Polya number is: 9784472371200 = 6!^2 * 4!^2 * 2!^15
The 2800th Jordan-Polya number is: 439378587648000 = 14! * 7!
The 3800th Jordan-Polya number is: 7213895789838336 = 4!^8 * 2!^16

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

function aupto(limit::T) where T <: Integer
    res = map(factorial, T(1):T(18))
    k = 2
    while k < length(res)
        rk = res[k]
        for j = 2:length(res)
            kl = res[j] * rk
            kl > limit && break
            while kl <= limit && kl  res
                push!(res, kl)
                kl *= rk
             end
        end
        k += 1
    end
    return sort!((sizeof(T) > sizeof(Int) ? T : Int).(res))[begin+1:end]
end

const factorials = map(factorial, 2:18)

""" Factor a J-P number into a smallest vector of factorials and their powers """
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"]
""" Express a positive integer as Unicode superscript digit characters """
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])!$(t[2] > 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!¹⁶

Kotlin

Translation of: Java
import java.util.*

object JordanPolyaNumbers {
    private val jordanPolyaSet = TreeSet<Long>()
    private val decompositions = HashMap<Long, TreeMap<Int, Int>>()

    @JvmStatic
    fun main(aArgs: Array<String>) {
        createJordanPolya()

        val belowHundredMillion = jordanPolyaSet.floor(100_000_000L)
        val jordanPolya = ArrayList(jordanPolyaSet)

        println("The first 50 Jordan-Polya numbers:")
        for (i in 0 until 50) {
            print(String.format("%5s%s", jordanPolya[i], if (i % 10 == 9) "\n" else ""))
        }
        println()

        println("The largest Jordan-Polya number less than 100 million: $belowHundredMillion")
        println()

        for (i in listOf(800, 1050, 1800, 2800, 3800)) {
            println("The $i th Jordan-Polya number is: ${jordanPolya[i - 1]} = ${toString(decompositions[jordanPolya[i - 1]]!!)}")
        }
    }

    private fun createJordanPolya() {
        jordanPolyaSet.add(1L)
        val nextSet = TreeSet<Long>()
        decompositions[1L] = TreeMap()
        var factorial = 1L

        for (multiplier in 2..20) {
            factorial *= multiplier
            val iterator = jordanPolyaSet.iterator()
            while (iterator.hasNext()) {
                var number = iterator.next()
                while (number <= Long.MAX_VALUE / factorial) {
                    val original = number
                    number *= factorial
                    nextSet.add(number)
                    decompositions[number] = TreeMap(decompositions[original]!!)
                    decompositions[number]?.merge(multiplier, 1) { a, b -> a + b }
                }
            }
            jordanPolyaSet.addAll(nextSet)
            nextSet.clear()
        }
    }

    private fun toString(aMap: Map<Int, Int>): String {
        return aMap.entries.joinToString(separator = " * ") { (key, value) ->
            "$key!${if (value == 1) "" else "^$value"}"
        }
    }
}

fun main(args: Array<String>) {
    JordanPolyaNumbers.main(arrayOf<String>())
}
Output:
The 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 less than 100 million: 99532800

The 800 th Jordan-Polya number is: 18345885696 = 2!^2 * 4!^7
The 1050 th Jordan-Polya number is: 139345920000 = 2! * 5!^3 * 8!
The 1800 th Jordan-Polya number is: 9784472371200 = 2!^15 * 4!^2 * 6!^2
The 2800 th Jordan-Polya number is: 439378587648000 = 7! * 14!
The 3800 th Jordan-Polya number is: 7213895789838336 = 2!^16 * 4!^8

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!)¹⁶

Mathematica / Wolfram Language

Bonus task TBD.

ClearAll[JPNumbers, largestUnder100M];

JPNumbers [n_Integer?Positive] := 
  Module[{factorials = Table[k!, {k, 1, n}]},
   Return [FixedPoint[Select[Union @@ Outer[Times, factorials, #], (# <= n! &)] &, {1}]];
];

jpNumbers = JPNumbers[19];
largestUnder100M = Max@Select[# < 10^8 &][jpNumbers];

Print["First 50 Jordan-Pólya numbers:"];
Print[StringJoin@Riffle[Partition[StringPadLeft[#, 8] & /@ ToString /@ jpNumbers[[1 ;; 50]], UpTo[10]], "\n"]];
Print["The largest Jordan-Pólya number less than 100 million: ", largestUnder100M];
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-Pólya number less than 100 million: 99532800

Pascal

Free Pascal

succesive add of next factorial in its power.keep sorted and without doublettes.
Now using Uint64 and only marking which factorial is used, which is unnecessary, but why not. It makes output easier
Runtime for TIO.RUN 127 ms (too short to be significant ) @home 2 ms
Using alternative dblLimit := 1 shl 6 as starting Limit and increase it by a factor of 256
This gets if MaxIdx = 3800 in 4ms.

program Jordan_Polya_Num;
{$IFDEF FPC}{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$COPERATORS ON}{$ENDIF}
{$IFDEF Windows}{$APPTYPE CONSOLE}{$ENDIF}
uses
  sysutils;
const
  MaxIdx = 3800;//7279 < 2^62
  maxFac = 25;//21!> 2^63
type
  tnum = Uint64;
  tpow= set of 0..31;// 1==2!^? ,2=3!^? 3=2!^?*3!^?
  tFac_mul = packed record
               fm_num : tnum;
               fm_pow : tpow;
               fm_high_idx : word;
               fm_high_pow : word;
             end;
  tpFac_mul = ^tFac_mul;
  tFacMulPow = array of tFac_mul;
  tFactorial = array[0..maxFac-2] of tnum;

var
  FacMulPowGes : tFacMulPow;
  Factorial: tFactorial;
  LastSearchFor :tFac_mul;
  dblLimit : tnum;

function CommatizeUint64(num:Uint64):AnsiString;
var
   fromIdx,toIdx :Int32;
Begin
  str(num,result);
  fromIdx := length(result);
  toIdx := fromIdx-1;
  if toIdx < 3 then
    exit;

  toIdx := 4*(toIdx DIV 3)+toIdx MOD 3 +1 ;
  setlength(result,toIdx);
  repeat
    result[toIdx]   := result[FromIdx];
    result[toIdx-1] := result[FromIdx-1];
    result[toIdx-2] := result[FromIdx-2];
    result[toIdx-3] := ',';
    dec(toIdx,4);
    dec(FromIdx,3);
  until FromIdx<=3;
end;

procedure Out_MulFac(idx:Uint32;const fm:tFac_mul);
var
  fac,
  num : tNum;
  FacIdx,pow : integer;
begin
  num := fm.fm_num;
  FacIdx := fm.fm_high_idx;
  write(CommatizeUint64(num):25,' = ');

  repeat
    pow := 0;
    fac := Factorial[FacIdx];
    while (num>=fac) AND (num mod Fac = 0) do
    Begin
      num := num DIV Fac;
      inc(pow);
    end;
    if pow = 0 then
      write(' 1')
    else
      if pow = 1 then
        write(' ',FacIdx+2,'!')
      else
        write(' (',FacIdx+2,'!)^',pow);
    if num = 1 then
      BREAK;
    repeat
      dec(FacIdx);
    until(FacIdx<0) OR (FacIdx in fm.fm_pow);
  until FacIdx < 0;
  writeln;

end;

procedure Out_I_th(i: integer);
begin
  if i < 0 then
    write(i:8,' too small');
  if i <= High(FacMulPowGes) then
  begin
    write(i:6,'-th : ');
    Out_MulFac(i,FacMulPowGes[i-1])
  end
  else
    writeln('Too big');
end;

procedure Out_First_N(n: integer);
var
  s,fmt : AnsiString;
  i,tmp : integer;
Begin
  if n<1 then
    EXIT;
  writeln('The first ',n,' Jordan-Polia numbers');
  s := '';
  If n > Length(FacMulPowGes) then
    n := Length(FacMulPowGes);
  dec(n);
  tmp := length(CommatizeUint64(FacMulPowGes[n].fm_num))+1;
  fmt := '%'+IntToStr(tmp)+'s';
  tmp := 72 DIV tmp;
  For i := 0 to n do
  Begin
    s += Format(fmt,[CommatizeUint64(FacMulPowGes[i].fm_num)]);
    if (i+1) mod tmp = 0 then
    Begin
      writeln(s);
      s := '';
    end;
  end;
  if s <>'' then
    writeln(s);
  writeln;
end;

procedure Initfirst;
var
  fac: tnum;
  i,j,idx: integer;
Begin
  fac:= 1;
  j := 1;
  idx := 0;
  For i := 2 to maxFac do
  Begin
    repeat
      inc(j);
      fac *= j;
    until j = i;
    Factorial[idx] := fac;
    inc(idx);
  end;
  Fillchar(LastSearchFor,SizeOf(LastSearchFor),#0);
  LastSearchFor.FM_NUM := 0;
//  dblLimit := 1 shl 53;
  dblLimit := 1 shl 5;
end;

procedure ResetSearch;
Begin
  setlength(FacMulPowGes,0);
end;

procedure GenerateFirst(idx:NativeInt;var res:tFacMulPow);
//generating the first entry with (2!)^n
var
  Fac_mul :tFac_mul;
  facPow,Fac : tnum;
  i,MaxPowOfFac : integer;
begin
  fac := Factorial[idx];
  MaxPowOfFac := trunc(ln(dblLimit)/ln(Fac))+1;
  setlength(res,MaxPowOfFac);

  with Fac_Mul do
  begin
    fm_num := 1;
    fm_pow := [0];
    fm_high_idx := 0;
  end;

  res[0] := Fac_Mul;
  facPow := 1;
  i := 1;
  repeat
    facPow *= Fac;
    if facPow >dblLimit then
      BREAK;
    with Fac_Mul do
    begin
      fm_num := facPow;
      fm_high_pow := i;
    end;
    res[i] := Fac_Mul;
    inc(i);
  until i = MaxPowOfFac;
  setlength(res,i);
end;

procedure DelDoublettes(var FMP:tFacMulPow);
//throw out doublettes,
var
  pNext,pCurrent : tpFac_mul;
  i, len,idx : integer;
begin
  len := 0;
  pCurrent := @FMP[0];
  pNext := pCurrent;
  For i := 0 to High(FMP)-1 do
  begin
    inc(pNext);
    // don't increment pCurrent if equal
    // pCurrent gets or stays the highest Value in n!^high_pow
    if pCurrent^.fm_num = pNext^.fm_num then
    Begin
      idx := pCurrent^.fm_high_idx;
      if idx < pNext^.fm_high_idx then
        pCurrent^  := pNext^
      else
        if idx = pNext^.fm_high_idx then
          if pCurrent^.fm_high_pow < pNext^.fm_high_pow then
            pCurrent^  := pNext^;
    end
    else
    begin
      inc(len);
      inc(pCurrent);
      pCurrent^  := pNext^;
    end;
  end;
  setlength(FMP,len);
end;

procedure QuickSort(var AI: tFacMulPow; ALo, AHi: Int32);
var
  Tmp :tFac_mul;
  Pivot : tnum;
  Lo, Hi : Int32;
begin
  Lo := ALo;
  Hi := AHi;
  Pivot := AI[(Lo + Hi) div 2].fm_num;
  repeat
    while AI[Lo].fm_num < Pivot do
      Inc(Lo);
    while AI[Hi].fm_num > Pivot do
      Dec(Hi);
    if Lo <= Hi then
    begin
      Tmp := AI[Lo];
      AI[Lo] := AI[Hi];
      AI[Hi] := Tmp;
      Inc(Lo);
      Dec(Hi);
    end;
  until Lo > Hi;
  if Hi > ALo then
    QuickSort(AI, ALo, Hi) ;
  if Lo < AHi then
    QuickSort(AI, Lo, AHi) ;
end;

function InsertFacMulPow(var res:tFacMulPow;Facidx:integer):boolean;
var
  Fac,FacPow,NewNum,limit : tnum;
  l_res,l_NewMaxPow,idx,i,j : Integer;
begin
  fac := Factorial[Facidx];
  if fac>dblLimit then
    EXIT(false);

  if length(res)> 0 then
  begin
    l_NewMaxPow := trunc(ln(dblLimit)/ln(Fac))+1;
    l_res := length(res);
    //calc new length, reduces allocation of big memory chunks
    //first original length + length of the new to insert
    j := l_res+l_NewMaxPow;
    //find the maximal needed elements which stay below  dbllimit
    // for every Fac^i
    idx := High(res);
    FacPow := Fac;
    For i := 1 to l_NewMaxPow do
    Begin
      limit := dblLimit DIV FacPow;
      if limit < 1 then
        BREAK;
      //search for the right position
      repeat
        dec(idx);
      until res[idx].fm_num<=limit;
      inc(j,idx);
      FacPow *= fac;
    end;
    j += 2;
    setlength(res,j);

    idx := l_res;
    FacPow := fac;
    For j := 1 to l_NewMaxPow do
    begin
      For i := 0 to l_res do
      begin
        NewNum := res[i].fm_num*FacPow;
        if NewNum>dblLimit then
          Break;
        res[idx]:= res[i];
        with res[idx] do
        Begin
          fm_num := NewNum;
          include(fm_pow,Facidx);
          fm_high_idx := Facidx;
          fm_high_pow := j;
        end;
        inc(idx);
      end;
      FacPow *= fac;
    end;
    setlength(res,idx);
    QuickSort(res,Low(res),High(res));
    DelDoublettes(res);
  end
  else
    GenerateFirst(Facidx,res);
  Exit(true);
end;

var
  i : integer;
BEGIn
  InitFirst;

  repeat
    ResetSearch;
    i := 0;
    repeat
      if Not(InsertFacMulPow(FacMulPowGes,i)) then
        BREAK;
      inc(i);
    until i > High(Factorial);
    //check if MaxIdx is found
    if (Length(FacMulPowGes) > MaxIdx) then
    begin
      if (LastSearchFor.fm_num<> FacMulPowGes[MaxIdx-1].fm_num) then
      Begin
        LastSearchFor := FacMulPowGes[MaxIdx-1];
        //the next factorial is to big, so search is done
        if LastSearchFor.fm_num < Factorial[i] then
          break;
      end
      else
        Break;
    end;
    if dblLimit> HIGH(tNUm) DIV 256 then
      BREAK;
    dblLimit *= 256;
  until false;

  write('Found ',length(FacMulPowGes),' Jordan-Polia numbers ');
  writeln('up to ',CommatizeUint64(dblLimit));
  writeln;

  Out_First_N(50);

  write('The last < 1E8 ');
  for i := 0 to High(FacMulPowGes) do
    if FacMulPowGes[i].fm_num > 1E8 then
    begin
      Out_MulFac(i,FacMulPowGes[i-1]);
      BREAK;
    end;
  writeln;

  Out_I_th(1);
  Out_I_th(100);
  Out_I_th(800);
  Out_I_th(1050);
  Out_I_th(1800);
  Out_I_th(2800);
  Out_I_th(3800);
END.
@home:
Found 3876 Jordan-Polia numbers up to 9,007,199,254,740,992

The first 50 Jordan-Polia 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 1,024 1,152 1,296 1,440 1,536
 1,728 1,920 2,048 2,304 2,592 2,880 3,072 3,456 3,840 4,096 4,320 4,608
 5,040 5,184

The last < 1E8                99,532,800 =  (6!)^2 4! (2!)^3

     1-th :                         1 =  1
   100-th :                    92,160 =  6! (2!)^7
   800-th :            18,345,885,696 =  (4!)^7 (2!)^2
  1050-th :           139,345,920,000 =  8! (5!)^3 2!
  1800-th :         9,784,472,371,200 =  (6!)^2 (4!)^2 (2!)^15
  2800-th :       439,378,587,648,000 =  14! 7!
  3800-th :     7,213,895,789,838,336 =  (4!)^8 (2!)^16
real	0m0,004s user	0m0,004s sys	0m0,000s

Perl

Translation of: Raku
Library: ntheory
use strict;
use warnings;
use feature 'say';

use ntheory 'factorial';
use List::AllUtils <max firstidx>;

sub table { my $t = 10 * (my $c = 1 + length max @_); ( sprintf( ('%'.$c.'d')x@_, @_) ) =~ s/.{1,$t}\K/\n/gr }

sub Jordan_Polya {
    my $limit = shift;
    my($k,@JP) = (2);
    push @JP, factorial $_ for 0..18;

    while ($k < @JP) {
        my $rk = $JP[$k];
        for my $l (2 .. @JP) {
            my $kl = $JP[$l] * $rk;
            last if $kl > $limit;
            LOOP: {
                my $p = firstidx { $_ >= $kl } @JP;
                if    ($p  < $#JP and $JP[$p] != $kl) { splice @JP, $p, 0, $kl }
                elsif ($p == $#JP                   ) {   push @JP,        $kl }
                $kl > $limit/$rk ? last LOOP : ($kl *= $rk)
            }
        }
        $k++
    }
    shift @JP; return @JP
}

my @JP = Jordan_Polya 2**27;
say "First 50 Jordan-Pólya numbers:\n" . table @JP[0..49];
say 'The largest Jordan-Pólya number before 100 million: ' . $JP[-1 + firstidx { $_ > 1e8 } @JP];
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

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)
    --
    -- Subtract prime powers of factorials off the prime powers of the jp number,
    -- only for factorials that have the same high prime factor as the remainder,
    -- and only putting things back on the todo list if still viable. 
    -- Somewhat slowish, but at least it /is/ very thorough.
    --
    sequence p = prime_powers(jp)
    integer lp = length(p), 
            mp = p[lp][1],  -- (max prime factor)
            hf = get_prime(lp+1)-1 -- (high factorial)
    assert(mp = get_prime(lp))
    sequence ap = get_primes_le(mp), -- (all primes)
             fs = apply(tagset(hf),factorial),
             fp = apply(fs,prime_powers),
             pf = repeat(0,hf), -- (powers of factorials)
             todo = {{p,pf}},
             seen = {},
             result = {}
    while length(todo) do
        {{p,pf},todo} = {todo[1],todo[2..$]}
        for fdx,fpi in fp from 2 do
            if fpi[$][1] = p[$][1] then -- same max prime factor
                bool ok = true
                for j,fpij in fpi do
                    if fpij[2]>p[find(fpij[1],ap)][2] then
                        -- this factorial ain't a factor
                        ok = false
                        exit
                    end if
                end for
                if ok then
                    -- reduce & trim the remaining prime powers:
                    sequence pnxt = deep_copy(p)
                    for j,fpij in fpi do
                        pnxt[find(fpij[1],ap)][2] -= fpij[2]
                    end for
                    while length(pnxt) and pnxt[$][2]=0 do
                        pnxt = pnxt[1..$-1]
                    end while
                    sequence fnxt = deep_copy(pf)
                    fnxt[fdx] += 1 -- **one** extra factorial power
                    if length(pnxt) then
                        bool bBad = false
                        for i=2 to length(pnxt) do
                            if pnxt[i][2]>pnxt[i-1][2] then
                                -- ie/eg you cannot ever knock a 7! or above off
                                --  if there ain't enough 5 (and 3 and 2) avail.
                                bBad = true
                                exit
                            end if
                        end for
                        if not bBad
                        and not find({pnxt,fnxt},seen) then
                            seen = append(seen,{pnxt,fnxt})
                            todo = append(todo,{pnxt,fnxt})
                        end if
                    else
                        result = append(result,fnxt)
                    end if
                end if
            end if
        end for
    end while       
    result = reverse(sort(apply(result,reverse))[$])
    string res = ""
    for i=length(result) to 1 by -1 do
        if result[i] then
            if length(res) then res &= " * " end if
            res &= sprintf("%d!",i)
            if result[i]>1 then
                res &= sprintf("^%d",result[i])
            end if
        end if
    end for
    return res
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,1050,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 1050th is 139,345,920,000 = 8! * 5!^3 * 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
"1.5s"

Some 80%-90% of the time is now spent in the decomposing phase.

Python

Translation of: Java
from collections import defaultdict
from itertools import product

class JordanPolyaNumbers:
    def __init__(self):
        self.jordan_polya_set = set()
        self.decompositions = defaultdict(dict)

    def create_jordan_polya(self):
        self.jordan_polya_set.add(1)
        next_set = set()
        self.decompositions[1] = {}
        factorial = 1

        for multiplier in range(2, 21):
            factorial *= multiplier
            for number in list(self.jordan_polya_set):
                while number <= 2**63 - 1 // factorial:
                    original = number
                    number *= factorial
                    next_set.add(number)
                    self.decompositions[number] = self.decompositions[original].copy()
                    self.decompositions[number][multiplier] = self.decompositions[number].get(multiplier, 0) + 1

            self.jordan_polya_set.update(next_set)
            next_set.clear()

    def to_string(self, a_map):
        result = ""
        for key in sorted(a_map.keys(), reverse=True):
            exponent = a_map[key]
            result += f"{key}!" + ("" if exponent == 1 else f"^{exponent}") + " * "
        return result[:-3]

    def display_results(self):
        below_hundred_million = max(n for n in self.jordan_polya_set if n < 100_000_000)
        jordan_polya = sorted(list(self.jordan_polya_set))

        print("The first 50 Jordan-Polya numbers:")
        for i in range(50):
            end = "\n" if (i % 10 == 9) else ""
            print(f"{jordan_polya[i]:5}", end=end)
        print()

        print(f"The largest Jordan-Polya number less than 100 million: {below_hundred_million}")
        print()

        for i in [800, 1050, 1800, 2800, 3800]:
            print(f"The {i}th Jordan-Polya number is: {jordan_polya[i-1]}"
                  f" = {self.to_string(self.decompositions[jordan_polya[i-1]])}")


if __name__ == "__main__":
    jpn = JordanPolyaNumbers()
    jpn.create_jordan_polya()
    jpn.display_results()
Output:
The 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 less than 100 million: 99532800

The 800th Jordan-Polya number is: 18345885696 = 4!^7 * 2!^2
The 1050th Jordan-Polya number is: 139345920000 = 8! * 5!^3 * 2!
The 1800th Jordan-Polya number is: 9784472371200 = 6!^2 * 4!^2 * 2!^15
The 2800th Jordan-Polya number is: 439378587648000 = 14! * 7!
The 3800th Jordan-Polya number is: 7213895789838336 = 4!^8 * 2!^16

Raku

Partial translation of Go

# 20230719 Raku programming solution

my \factorials = 1, | [\*] 1..18; # with 0!

sub JordanPolya (\limit) {
   my \ix = (factorials.keys.first: factorials[*] >= limit) // factorials.end; 
   my ($k, @res) = 2, |factorials[0..ix];

   while $k < @res.elems {
      my \rk = @res[$k];
      for 2 .. @res.elems -> \l {
         my \kl = $ = @res[l] * rk;
         last if kl > limit;
         loop {
            my \p = @res.keys.first: { @res[$_] >= kl } # performance
            if p < @res.elems and @res[p] != kl {
               @res.splice: p, 0, kl
            } elsif p == @res.elems { 
               @res.append: kl 
            }
            kl > limit/rk ?? ( last ) !! kl *= rk
         }
      }
      $k++
   }
   return @res[1..*]
}

my @result = JordanPolya 2**30 ; 
say "First 50 Jordan-Pólya numbers:";
say [~] $_>>.fmt('%5s') for @result[^50].rotor(10);
print "\nThe largest Jordan-Pólya number before 100 million: ";
say @result.first: * < 100_000_000, :end;

You may Attempt This Online!

Rust

Translation of: C++
use std::collections::{HashMap, HashSet};

const LIMIT: u64 = 1 << 53;

fn to_string(map: &HashMap<u64, u64>) -> String {
    let mut result = String::new();
    for (&key, &value) in map {
        // println!("key={} value={}", key, value);
        let part = if value == 1 {
            format!("{}!", key)
        } else {
            format!("{}!^{}", key, value)
        };
        result = part + " * " + &result;
    }
    result.trim_end_matches(" * ").to_string()
}

fn insert_or_update(map: &mut HashMap<u64, u64>, entry: u64) {
    let count = map.entry(entry).or_insert(0);
    *count += 1;
}

fn create_jordan_polya() -> (HashSet<u64>, HashMap<u64, HashMap<u64, u64>>) {
    let mut jordan_polya_set = HashSet::new();
    let mut decompositions = HashMap::new();

    jordan_polya_set.insert(1);
    decompositions.insert(1, HashMap::new());
    let mut factorial = 1u64;

    for multiplier in 2..=20 {
        factorial *= multiplier; // Using u64 for multiplier
        let mut to_update = Vec::new();

        for &number in &jordan_polya_set {
            let mut current = number;
            while current <= LIMIT / factorial {
                to_update.push((current, current * factorial)); // Store original and new number
                current *= factorial;
            }
        }

        for (original, new_number) in to_update {
            jordan_polya_set.insert(new_number);
            let mut new_decomposition = decompositions[&original].clone();
            insert_or_update(&mut new_decomposition, multiplier);
            decompositions.insert(new_number, new_decomposition);
        }
    }

    (jordan_polya_set, decompositions)
}


fn main() {
    let (jordan_polya_set, decompositions) = create_jordan_polya();
    let mut jordan_polya: Vec<_> = jordan_polya_set.into_iter().collect();
    jordan_polya.sort();

    println!("The first 50 Jordan-Polya numbers:");
    for i in 0..50 {
        print!("{:5}", jordan_polya[i]);
        if i % 10 == 9 {
            println!();
        }
    }

    let hundred_million = jordan_polya.iter().position(|&x| x >= 100_000_000).unwrap();
    println!(
        "\nThe largest Jordan-Polya number less than 100 million: {}\n",
        jordan_polya[hundred_million - 1]
    );

    for &i in &[800, 1050, 1800, 2800, 3800] {
        println!(
            "The {}th Jordan-Polya number is: {} = {}",
            i,
            jordan_polya[i - 1],
            to_string(&decompositions[&jordan_polya[i - 1]])
        );
    }
}
Output:
The 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 less than 100 million: 99532800

The 800th Jordan-Polya number is: 18345885696 = 2!^2 * 4!^7
The 1050th Jordan-Polya number is: 139345920000 = 2! * 5!^3 * 8!
The 1800th Jordan-Polya number is: 9784472371200 = 2!^15 * 6!^2 * 4!^2
The 2800th Jordan-Polya number is: 439378587648000 = 14! * 7!
The 3800th Jordan-Polya number is: 7213895789838336 = 2!^16 * 4!^8

Scala

Translation of: Java
import java.util.{ArrayList, HashMap, TreeMap, TreeSet}
import scala.jdk.CollectionConverters._

object JordanPolyaNumbers {
  
  private val jordanPolyaSet = new TreeSet[Long]()
  private val decompositions = new HashMap[Long, TreeMap[Integer, Integer]]()

  def main(args: Array[String]): Unit = {
    createJordanPolya()

    val belowHundredMillion = jordanPolyaSet.floor(100_000_000L)
    val jordanPolya = new ArrayList[Long](jordanPolyaSet)

    println("The first 50 Jordan-Polya numbers:")
    jordanPolya.asScala.take(50).zipWithIndex.foreach { case (number, index) =>
      print(f"$number%5s${if(index % 10 == 9) "\n" else ""}")
    }
    println()

    println(s"The largest Jordan-Polya number less than 100 million: $belowHundredMillion")
    println()

    List(800, 1050, 1800, 2800, 3800).foreach { i =>
      println(s"The ${i}th Jordan-Polya number is: ${jordanPolya.get(i - 1)} = ${toString(decompositions.get(jordanPolya.get(i - 1)))}")
    }
  }

  private def createJordanPolya(): Unit = {
    jordanPolyaSet.add(1L)
    val nextSet = new TreeSet[Long]()
    decompositions.put(1L, new TreeMap[Integer, Integer]())
    var factorial = 1L

    for (multiplier <- 2 to 20) {
      factorial *= multiplier
      val it = jordanPolyaSet.iterator()
      while (it.hasNext) {
        var number = it.next()
        while (number <= Long.MaxValue / factorial) {
          val original = number
          number *= factorial
          nextSet.add(number)
          decompositions.put(number, new TreeMap[Integer, Integer](decompositions.get(original)))
          val currentMap = decompositions.get(number)
          currentMap.merge(multiplier, 1, (a: Integer, b: Integer) => Integer.sum(a, b))
        }
      }
      jordanPolyaSet.addAll(nextSet)
      nextSet.clear()
    }
  }

  private def toString(aMap: TreeMap[Integer, Integer]): String = {
    aMap.descendingMap().asScala.map { case (key, value) =>
      s"$key!${if (value == 1) "" else "^" + value}"
    }.mkString(" * ")
  }
}
Output:
The 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 less than 100 million: 99532800

The 800th Jordan-Polya number is: 18345885696 = 4!^7 * 2!^2
The 1050th Jordan-Polya number is: 139345920000 = 8! * 5!^3 * 2!
The 1800th Jordan-Polya number is: 9784472371200 = 6!^2 * 4!^2 * 2!^15
The 2800th Jordan-Polya number is: 439378587648000 = 14! * 7!
The 3800th Jordan-Polya number is: 7213895789838336 = 4!^8 * 2!^16

You may Attempt This Online!

Swift

Translation of: Java
import Foundation

class JordanPolyaNumbers {

    static var jordanPolyaSet = Set<Int64>()
    static var decompositions = [Int64: [Int: Int]]()

    static func main() {
        createJordanPolya()

        let belowHundredMillion = jordanPolyaSet.filter { $0 <= 100_000_000 }.max() ?? 0
        let jordanPolya = Array(jordanPolyaSet).sorted()

        print("The first 50 Jordan-Polya numbers:")
        for i in 0..<50 {
            if i % 10 == 9 {
                print("\(jordanPolya[i])\n", terminator: "")
            } else {
                print(String(format: "%5d ", jordanPolya[i]), terminator: "")
            }
        }
        print("\nThe largest Jordan-Polya number less than 100 million: \(belowHundredMillion)\n")

        for i in [800, 1050, 1800, 2800, 3800] {
            print("The \(i)th Jordan-Polya number is: \(jordanPolya[i - 1]) = \(toString(decompositions[jordanPolya[i - 1]] ?? [:]))")
        }
    }

    static func createJordanPolya() {
        jordanPolyaSet.insert(1)
        var nextSet = Set<Int64>()
        decompositions[1] = [:]
        var factorial: Int64 = 1

        for multiplier in 2...20 {
            factorial *= Int64(multiplier)
            for number in jordanPolyaSet {
                var current = number
                while current <= Int64.max / factorial {
                    let original = current
                    current *= factorial
                    nextSet.insert(current)
                    decompositions[current] = decompositions[original]?.merging([multiplier: 1], uniquingKeysWith: +) ?? [:]
                }
            }
            jordanPolyaSet.formUnion(nextSet)
            nextSet.removeAll()
        }
    }

    static func toString(_ aMap: [Int: Int]) -> String {
        var result = ""
        for key in aMap.keys.sorted().reversed() {
            let value = aMap[key] ?? 0
            result += "\(key)!" + (value == 1 ? "" : "^\(value)") + " * "
        }
        return String(result.dropLast(3))
    }
}

JordanPolyaNumbers.main()
Output:
The 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 less than 100 million: 99532800

The 800th Jordan-Polya number is: 18345885696 = 4!^7 * 2!^2
The 1050th Jordan-Polya number is: 139345920000 = 8! * 5!^3 * 2!
The 1800th Jordan-Polya number is: 9784472371200 = 6!^2 * 4!^2 * 2!^15
The 2800th Jordan-Polya number is: 439378587648000 = 14! * 7!
The 3800th Jordan-Polya number is: 7213895789838336 = 4!^8 * 2!^16

Wren

Version 1

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 = []
    while (m % factorials[start] == 0) {
        f.add(start)
        m =  m / factorials[start]
        if (m == 1) return f
    }
    if (f.count > 0) {
        var g = Decompose.call(m, start-1)
        if (g.count > 0) {
            var prod = 1
            for (e in g) prod = prod * factorials[e]
            if (prod == m) return f + g
        }
    }
    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, 1050, 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,050th Jordan-Pólya number is : 139,345,920,000
= 8! x (5!)³ 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!)¹⁶

Version 2

Library: Wren-sort

This uses the same non-recursive algorithm as the Phix entry to generate the J-P numbers which, at 1.1 seconds on my machine, is about 40 times quicker than the OEIS algorithm.

import "./sort" for Find
import "./seq" for Lst
import "./fmt" for Fmt

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 JordanPolya = Fn.new { |limit|
    var ix = Find.nearest(factorials, limit).min(18)
    var res = factorials[0..ix]
    var k = 2
    while (k < res.count) {
        var rk = res[k]
        for (l in 2...res.count) {
            var kl = res[l] * rk
            if (kl > limit) break
            while (true) {
                var p = Find.nearest(res, kl)
                if (p < res.count && res[p] != kl) {
                    res.insert(p, kl)
                } else if (p == res.count) {
                    res.add(kl)
                }
                kl = kl * rk
                if (kl > limit) break
            }
        }
        k = k + 1
    }
    return res[1..-1]
}

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

cacheFactorials.call()
var v = JordanPolya.call(2.pow(53)-1)
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, 1050, 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:
Identical to first version.

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