How to iteratively build a list?

Question

I have the following c++ code that I want to translate to Mathematica.

std::vector<int> squbes;
for (int p : primes) for (int q : primes) if (p != q){
    if (p*p*q*q*q > maxi) break;
    squbes.push_back(p*p*q*q*q);
}

As you can see I simply have a list of numbers that I want to iteratively append elements to.

The first thing I tried was the following functional approach.

sqube[{x_,y_}] := x*x*y*y*y
unequal[{x_, y_}] := x != y
smallenough[x_] := x <= maxi
pairs := Select[Tuples[primes, 2], unequal]
squbes := Select[Map[sqube, pairs], smallenough]

However this is way too slow because it misses the very important break condition in the loop in the c++ version.

Next I tried the most direct translation I could come up with.

sqube[{x_,y_}] := x*x*y*y*y
squbes := {}
For[i = 1, i <= Length[primes], i++,
    For[j = 1, j <= Length[primes] && sqube[{primes[[i]], primes[[j]]}] <= maxi, j++, 
        AppendTo[squbes, sqube[{primes[[i]], primes[[j]]}]]]]

However for some reason this turned out to be really slow, even if primes has only $100$ elements it took more than a second. I suspect AppendTo creates a completely new list every time?

The final thing I tried was to use ReplacePart, but the attempt is not even worth showing. First problem is that you need to initialize squbes to be large enough, even though I have no good idea how large it will be. Second problem is that I see no reason to assume that ReplacePart will not create a completely new list every time.

Can you please help me simply create this list of squbes in a decent amount of time?

Answer 1

Let me see if I get this right. So you have

max = 10^12

and you are iterating over the same list of primes with a nested loop. Your inner loop iterates over q and we should first calculate the largest prime we need in the list to ensure that we are able to hit max.

Your very first iteration, where p=2 is critical because p^2*q^3 will be the smallest. So how about we find out for which prime q we make it above the max mark? Just for clarity, I'm making the code a bit verbose

With[{p = 2},
 MinimalBy[Range[2000],
  Abs[p^2*Prime[#]^3 - max] &
  ]
 ]
(* {819} *)

This tells us that we are passing the max mark for the 820th prime

2^2*Prime[820]^3
(* 1000664355604 *)

For collecting the values, I suggest you use Reap and Sow instead of AppendTo as it is faster. For the iteration, there are many choices. Let me use a simple Do in this case

result = Reap[With[
  {
    primes = Prime[Range[820]]
  },
  Do[
    If[p =!= q && p^2*q^3 < max,
      Sow[p^2*q^3]
    ],
    {p, primes},
    {q, primes}
  ]
]][[2, 1]];

That last [[2,1]] looks a bit clumsy, but when you read how Reap works, you'll understand its meaning. I got about 21k results and they all have the required form

result[[500]] // FactorInteger
(* {{2, 2}, {3581, 3}} *)

Answer 2

The GeneralUtilities`` package contains undocumented iterator functionality that resembles the generator-based iterators from C++. Perhaps these will make it into the core language some day. Until then, the following iterator-based answer is mainly out of academic interest rather than a practical recommendation...

We start by defining sqube as in the question:

sqube[{x_,y_}] := x*x*y*y*y

We then define an iterator over all the primes (RangeIterator does not support Infinity so we use a googol instead):

primes[] := RangeIterator[1*^100] // MapIterator[Prime]

So then:

primes[] // TakeIterator[10] // Normal

(* {2, 3, 5, 7, 11, 13, 17, 19, 23, 29} *)

We will also define a useful helper iterator (which really ought to be built-in):

takeWhileIterator[crit_] := MapIterator[If[crit[#], #, IteratorExhausted]&]

Armed with these definitions, we can define our results:

results[maxi_] :=
  (* for (int p : primes) *)
  primes[] //
  (* not present in the C++, removes the need to know in advance how many primes *)
  takeWhileIterator[sqube[{2, #}] <= maxi &] //
  (* for (int q : primes) *)
  JoinMapIterator[
     ( primes[] //
       (* if (p != q) *)
       SelectIterator[Curry[#2 != # &][#]] //
       (* compute p*p*q*q*q *)
       MapIterator[Curry[sqube[{#2, #}]&][#]] //
       (* if (p*p*q*q*q > maxi) break *)
       takeWhileIterator[# <= maxi &]
     ) &
  ]

So then:

results[500] // Normal
(* {108, 500, 72, 200} *)

results[5000] // Normal
(* {108, 500, 1372, 72, 1125, 3087, 200, 675, 392, 1323} *)

results[10^12] // Normal // Short
(* { 108, 500, 1372, 5324, 8788, 19652, 27436
   , <<20952>>
   , 52810620731, 87171249997, 194935071113
   , 272147293459, 482754937967, 967692132989
   }
*)

% // Last // FactorInteger
(* {{29, 3}, {6299, 2}} *)