Like I commented, every single number required by the OP has exactly six divisors (including one and itself, with the trivial exception of 1). Consequentially, there are only two possible options: either the number is a fifth power of a prime: $p^5$ and has divisors $p^{0..5}$ whose product is, of course, $p^{15} = {(p^5)}^3$; otherwise it is the product $p_i^2 p_j$ and has divisors $1, p_i, p_j, p_i p_j, p_i^2, p_i^2 p_j$. Their product is of course $p_i^6 p_j^3 = {(p_i^2 p_j)}^3$.
I now want to generate all possible tuples of the form {i, i, j}
for i<=n
and j<=n
(and i != j
). These kind of tasks are best met by built-in functions like Subsets
, but since the first number is repeated, a customized compiled function to suit my needs should also work very well.
Code dump:
getTuples =
Compile[{{n, _Integer}},
Block[{out = Table[0, {n * (n - 1)}, {3}], i = 1, j = 1, k = 1},
For[i = 1, i <= n, i++,
For[j = 1, j < i, j++; k++,
out[[k, {1, 2}]] = Table[i, {2}]; out[[k, 3]] = j];
For[j = i + 1, j <= n, j++; k++,
out[[k, {1, 2}]] = Table[i, {2}]; out[[k, 3]] = j]
]; out]]
to actually get the required numbers we run the following:
getTuples[10] // Prime // Times @@@ # &
We lack here the fifth powers of primes though. To remedy, I define a little helper function:
getNumbers[n_] :=
Sort[{1}~Join~((Prime[Range[n]])^5)~
Join~(getTuples[n] // Prime // Times @@@ # &)]
It gets me 1000001 of such numbers in 3 seconds, however there is no guarantee that they are the first 1000001; since $p_i^2 p_j$ for $j = 1, i = 1000$ is present, but the likely much smaller case of, say, $j = 1001, i = 2$ is not.
Obligatory sample output:
getNumbers[10]
(* {1, 12, 18, 20, 28, 32, 44, 45, 50, 52, 63, 68, 75, 76, 92, 98, 99, \
116, 117, 147, 153, 171, 175, 207, 242, 243, 245, 261, 275, 325, 338, \
363, 425, 475, 507, 539, 575, 578, 605, 637, 722, 725, 833, 845, 847, \
867, 931, 1058, 1083, 1127, 1183, 1421, 1445, 1573, 1587, 1682, 1805, \
1859, 2023, 2057, 2299, 2523, 2527, 2645, 2783, 2873, 3125, 3179, \
3211, 3509, 3703, 3757, 3887, 3971, 4205, 4693, 4901, 5491, 5819, \
5887, 6137, 6647, 6877, 8303, 8381, 8993, 9251, 10051, 10469, 10933, \
14297, 15341, 15979, 16807, 19343, 161051, 371293, 1419857, 2476099, \
6436343, 20511149} *)
Observe how it fails though:
ListPlot[{getNumbers[10], getNumbers[100][[;; 101]]}, PlotRange -> {Full, {0, 1000}}]

The output of the first 19 primes matches for both lists, but then getNumbers[10]
shoots up, missing several suitable values. There may be ways to improve this.