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.
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. 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
- Wikipedia article : Jordan-Pólya number
- OEIS sequence A001013: Jordan-Pólya numbers
C
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!)¹⁶
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
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!¹⁶
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!)¹⁶
Pascal
Free Pascal
succesive add of next factorial in its power.keep sorted and without doublettes.
I dont't know, how "far" extended gets correct results.Maybe using logs would be more precise.
program Jordan_Polya_Num;
{$IFDEF FPC}{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$ENDIF}
{$IFDEF Windows}{$APPTYPE CONSOLE}{$ENDIF}
uses
sysutils;
const
dblLimit = 1E19;//7213895789838336+3e14;//1e53;
maxFac = 43;
type
tnum = extended;
tpow= array[0..maxFac-2] of byte;
tFac_mul = packed record
fm_num : tnum;
fm_pow : tpow;
fm_idx : byte;
end;
tpFac_mul = ^tFac_mul;
tFacMulPow = array of tFac_mul;
var
Factorial: array[0..maxFac-2] of tnum;
FacMulPowGes : tFacMulPow;
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;
procedure Out_MulFac(const fm:tFac_mul);
var
i,j,pow : integer;
begin
if fm.fm_num < 1E20 then
write(fm.fm_num:20:0)
else
writeln(fm.fm_num);
i := High(tpow);
while (i>=0 ) AND (fm.fm_pow[i]= 0) do
dec(i);
For j := 0 to i do
Begin
pow := fm.fm_pow[j];
if pow > 1 then
write(' (',j+2,'!)^',pow)
else
if pow= 1 then
write(' ',j+2,'!');
end;
writeln;
end;
procedure Init;
var
fac: tnum;
i,j,idx: integer;
Begin
fac:= 1.0;
j := 1;
idx := 0;
For i := 2 to 43 do
Begin
repeat
inc(j);
fac *= j;
until j = i;
Factorial[idx] := fac;
inc(idx);
end;
end;
procedure GenerateFirst(idx:NativeInt;var res:tFacMulPow);
//generating the first entry with (2!)^n
var
Fac_mul :tFac_mul;
fac : tnum;
i,MaxIdx : integer;
begin
fac := Factorial[idx];
MaxIDx := trunc(ln(dblLimit)/ln(Fac))+1;
setlength(res,MaxIdx);
fillchar(Fac_Mul,SizeOf(Fac_Mul),#0);
with Fac_Mul do
begin
fm_num := 1;
fm_pow[idx] := 0;
fm_idx := 0;
end;
res[0] := Fac_Mul;
fac := 1;
For i := 1 to MaxIdx-1 do
begin
fac *= Factorial[idx];
with Fac_Mul do
begin
fm_num := fac;
fm_pow[idx] := i;
end;
res[i] := Fac_Mul;
end;
end;
procedure DelDoulettes(var FMP:tFacMulPow);
//throw out doublettes,
//the one with highest power in the highest n! survives
var
pI,pJ : tpFac_mul;
i, j,idx : integer;
begin
j := 0;
pJ := @FMP[0];
pI := pJ;
For i := 0 to High(FMP)-1 do
begin
inc(pI);
if pJ^.fm_num = pI^.fm_num then
Begin
idx := pJ^.fm_idx;
if idx < pI^.fm_idx then
pJ^ := pI^
else
if idx = pI^.fm_idx then
if pJ^.fm_pow[idx]<pI^.fm_pow[idx] then
pJ^ := pI^;
end
else
begin
inc(j);
inc(pJ);
pJ^ := pI^;
end;
end;
setlength(FMP,j);
end;
procedure InsertFacMulPow(var res:tFacMulPow;Facidx:integer);
var
Fac,newFac,limit : tnum;
l_res,l_NewMaxPow,idx,i,j : Integer;
begin
if length(res)= 0 then
Begin
GenerateFirst(Facidx,res);
EXIT;
end;
fac := Factorial[Facidx];
if fac>dblLimit then
EXIT;
l_NewMaxPow := trunc(ln(dblLimit)/ln(Fac))+1;
l_res := length(res);
//calc new length, reduces allocation of big memory chunks
j := 0;
idx := High(res);
For i := 1 to l_NewMaxPow do
Begin
limit := dblLimit/fac;
if limit < 1 then
BREAK;
repeat
dec(idx);
until res[idx].fm_num<=limit;
inc(j,idx);
fac *=Factorial[Facidx];
end;
j += l_res+l_NewMaxPow+2;
setlength(res,j);
fac := Factorial[Facidx];
idx := l_res;
For j := 0 to l_NewMaxPow-1 do
begin
For i := 0 to l_res-1 do
begin
res[idx]:= res[i];
NewFac := res[i].fm_num*Fac;
if NewFac>dblLimit then
Break;
res[idx].fm_num := NewFac;
res[idx].fm_pow[Facidx] := j+1;
res[idx].fm_idx := Facidx;
inc(idx);
end;
fac *= Factorial[Facidx];
end;
setlength(res,idx);
QuickSort(res,Low(res),High(res));
DelDoulettes(res);
end;
var
i : integer;
BEGIn
init;
For i := Low(Factorial) to High(Factorial) do
InsertFacMulPow(FacMulPowGes,i);
write('Found ',length( FacMulPowGes),' Jordan-Polia numbers ');
writeln('up to ',dblLimit);
writeln;
writeln('The first 50 Jordan-Polia numbers');
For i := 1 to 50 do
Begin
write(FacMulPowGes[i-1].fm_num:5:0);
if i mod 10 = 0 then
writeln;
end;
writeln;
writeln('The last < 1E8 ');
for i := 0 to High(FacMulPowGes) do
if FacMulPowGes[i].fm_num >= 1E8 then
begin
write('Index: ',i,' = ');
Out_MulFac(FacMulPowGes[i-1]);
BREAK;
end;
writeln;
writeln(' Index ');
i := 100;write(i:8,': ');Out_MulFac(FacMulPowGes[i-1]);
i := 800;write(i:8,': ');Out_MulFac(FacMulPowGes[i-1]);
i := 1050;write(i:8,': ');Out_MulFac(FacMulPowGes[i-1]);
i := 1800;write(i:8,': ');Out_MulFac(FacMulPowGes[i-1]);
i := 2800;write(i:8,': ');Out_MulFac(FacMulPowGes[i-1]);
i := 3800;write(i:8,': ');Out_MulFac(FacMulPowGes[i-1]);
END.
- @TIO.RUN:
Found 7832 Jordan-Polia numbers up to 1.00000000000000000000E+0019
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 1024 1152 1296 1440 1536 1728 1920 2048 2304
2592 2880 3072 3456 3840 4096 4320 4608 5040 5184
The last < 1E8
Index: 367 = 99532800 (2!)^3 4! (6!)^2
Index
100: 92160 (2!)^7 6!
800: 18345885696 (2!)^2 (4!)^7
1050: 139345920000 2! (5!)^3 8!
1800: 9784472371200 (2!)^15 (4!)^2 (6!)^2
2800: 439378587648000 7! 14!
3800: 7222041363087360 2! (3!)^11 (4!)^3 6!
Real time: 0.148 s User time: 0.122 s Sys. time: 0.024 s CPU share: 99.01 %
Found 1660536 Jordan-Polia numbers up to 9.99999999999999999971E+0052
--using double Found 1933972 Jordan-Polia numbers up to 9.99999999999999999971E+0052
Real time: 13.738 s User time: 12.925 s Sys. time: 0.715 s CPU share: 99.28 %
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.
Wren
Version 1
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
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