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;

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?

  • $\begingroup$ an aside: N is a reserved symbol. You should get an error message if you try to assign a value to it. $\endgroup$ – kglr Jul 21 at 18:55
  • $\begingroup$ Fair point, I changed it to maxi. $\endgroup$ – SmileyCraft Jul 21 at 18:56
  • $\begingroup$ what is the interesting range for maxi? $\endgroup$ – kglr Jul 21 at 19:01
  • $\begingroup$ maxi is about $10^{12}$ $\endgroup$ – SmileyCraft Jul 21 at 19:03
  • 2
    $\begingroup$ That's right. But he who appends things to a list in a long loop needs to have good reason for that, because it is allways a major hit to performance. If you really require dynamic appending, then I would suggest to get to know to (i) Association, (ii) Sow and Reap, or (iii) Internal`Bag. The latter is not documented, though, but you can find details about it on this site. $\endgroup$ – Henrik Schumacher Jul 21 at 19:28

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},
  Abs[p^2*Prime[#]^3 - max] &
(* {819} *)

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

(* 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]]
    If[p =!= q && p^2*q^3 < max,
    {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}} *)

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) *)
     ( 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}} *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.