(* first we compute the primordial numbers up to 2^powerLim *)
powerLim = 20; pr = Prime[Range[powerLim]]; len = Length[pr];
t = Flatten[Table[pr^PadRight[p, len], {n, len}, {p, IntegerPartitions[n]}], 1];
t2 =Table[Times @@ i, {i, t}]; pri = Sort[Select[t2, # <= pr[[1]]^len &]]; Length[pri]
(* then we find the numbers that produce a record number of solutions *)
mxSols = 0;
tAll = Table[
fact = FactorInteger[theNum]; expon = Transpose[fact][[2]];
len = Length[fact];
Switch[len,
1, t = Flatten[Table[2^a,
{a, Compositions[expon[[1]], 3]}], 0];,
2, t = Flatten[Table[2^a 3^b,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]}], 1];,
3, t = Flatten[Table[2^a 3^b 5^c,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]}], 2];,
4, t = Flatten[Table[2^a 3^b 5^c 7^d,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]}], 3];,
5, t = Flatten[Table[2^a 3^b 5^c 7^d 11^e,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]},
{e, Compositions[expon[[5]], 3]}], 4];,
6, t = Flatten[Table[2^a 3^b 5^c 7^d 11^e 13^f,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]},
{e, Compositions[expon[[5]], 3]},
{f, Compositions[expon[[6]], 3]}], 5];,
7, t = Flatten[Table[2^a 3^b 5^c 7^d 11^e 13^f 17^g,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]},
{e, Compositions[expon[[5]], 3]},
{f, Compositions[expon[[6]], 3]},
{g, Compositions[expon[[7]], 3]}], 6];,
8, t = Flatten[Table[2^a 3^b 5^c 7^d 11^e 13^f 17^g 19^h,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]},
{e, Compositions[expon[[5]], 3]},
{f, Compositions[expon[[6]], 3]},
{g, Compositions[expon[[7]], 3]},
{h, Compositions[expon[[8]], 3]}], 7];,
9, t = Flatten[Table[2^a 3^b 5^c 7^d 11^e 13^f 17^g 19^h 23^i,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]},
{e, Compositions[expon[[5]], 3]},
{f, Compositions[expon[[6]], 3]},
{g, Compositions[expon[[7]], 3]},
{h, Compositions[expon[[8]], 3]},
{i, Compositions[expon[[9]], 3]}], 8];,
10, t = Flatten[Table[2^a 3^b 5^c 7^d 11^e 13^f 17^g 19^h 23^i 29^j,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]},
{e, Compositions[expon[[5]], 3]},
{f, Compositions[expon[[6]], 3]},
{g, Compositions[expon[[7]], 3]},
{h, Compositions[expon[[8]], 3]},
{i, Compositions[expon[[9]], 3]},
{j, Compositions[expon[[10]], 3]}], 9];,
11, t = Flatten[Table[2^a 3^b 5^c 7^d 11^e 13^f 17^g 19^h 23^i 29^j 31^k,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]},
{e, Compositions[expon[[5]], 3]},
{f, Compositions[expon[[6]], 3]},
{g, Compositions[expon[[7]], 3]},
{h, Compositions[expon[[8]], 3]},
{i, Compositions[expon[[9]], 3]},
{j, Compositions[expon[[10]], 3]},
{k, Compositions[expon[[11]], 3]}], 10];,
12, t = Flatten[Table[2^a 3^b 5^c 7^d 11^e 13^f 17^g 19^h 23^i 29^j 31^k 37^l,
{a, Compositions[expon[[1]], 3]},
{b, Compositions[expon[[2]], 3]},
{c, Compositions[expon[[3]], 3]},
{d, Compositions[expon[[4]], 3]},
{e, Compositions[expon[[5]], 3]},
{f, Compositions[expon[[6]], 3]},
{g, Compositions[expon[[7]], 3]},
{h, Compositions[expon[[8]], 3]},
{i, Compositions[expon[[9]], 3]},
{j, Compositions[expon[[10]], 3]},
{k, Compositions[expon[[11]], 3]},
{l, Compositions[expon[[12]], 3]}], 11];];
t2 = Union[Sort /@ t];
t3 = Total /@ t2; t4 = Sort[t3];
t5 = Sort[Transpose[RotateLeft[Transpose[Tally[t4]], 1]]];
t6 = Reverse[Select[t5, PrimeQ[#[[2]]] &]];
If[Length[t6] == 0, {0, 0},
If[Length[t6] == 1,
If[t6[[1, 1]] > mxSols, mxSols = t6[[1, 1]];
Print[{theNum, t6[[1, 2]], expon, mxSols}]]; t6[[1]],
mx = t6[[1, 1]]; i = 2;
While[i <= Length[t6] && t6[[i, 1]] == mx, i++];
If[t6[[i - 1, 1]] > mxSols, mxSols = t6[[i - 1, 1]];
Print[{theNum, t6[[i - 1, 2]], expon, mxSols}]];
t6[[i - 1]]]], {theNum, pri}];