Can this be written well, without loops?

Question

Inspired by this question I would like to know if the following code can be written without explicit loops (For, While, etc.) in a clean, efficient and non-contrived way. I have been unable to do so.

max = 5000;
a = ConstantArray[0, max];
x = y = z = n = 1;
val := 2 (2 n^2+(y-2) (z-2)+x (y+z-2)+2 n (x+y+z-3));
For[x = 1, val <= max, x++,
 For[y = 1, val <= max && y <= x, y++,
  For[z = 1, val <= max && z <= y, z++,
   For[n = 1, (r = val) <= max, n++,
    a[[r]]++
   ]; n = 1
  ]; z = 1
 ]; y = 1
]

The output is the array a.

Answer 1

EDIT To address hard-coded Table and SparseArray limits, and efficiency

As pointed out in the comments, hard-coded limits on the Table or SparseArray dimensions may not work in general. Besides being slow, the Table approach quickly eats up system memory for moderate values of max. Here is a variation on WReach's recursive scheme using ReplaceRepeated. With max=5000, it is about a factor of 4 slower than using For.

Clear[max, a4];
max = 5000;
a4 = ConstantArray[0, max];
ReplaceRepeated[{1, 1, 1, 1},
 {
  {x_, y_, z_, n_} /; (r = 2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 2 n (x + y + z - 3))) 
                      <= max :> (If[z <= y <= x, a4[[r]]++]; {x, y, z, n + 1}),
  (* Stop *)
  {x_, 1, 1, 1} :> Null,
  (* Optimizations *)
  {x_, y_, 1, 1} :> If[y < x, {x, y + 1, 1, 1}, {x + 1, 1, 1, 1}],
  {x_, y_, z_, 1} :> If[z < y, {x, y, z + 1, 1}, {x, y + 1, 1, 1}],
  {x_, y_, z_, _} :> If[z < y, {x, y, z + 1, 1}, 
                        If[y < x, {x, y + 1, 1, 1}, {x + 1, 1, 1, 1}]]
 }
 , MaxIterations -> Infinity]

(Array-based solutions)

As far as readability, Table comes to mind:

Clear[val, a1, max];
max = 100;
a1 = ConstantArray[0, max];
val := 2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 2 n (x + y + z - 3));
Table[If[val <= max, a1[[val]]++], {x, 1, max}, {y, 1, x}, {z, 1, y}, {n, 1, max}];

a1==a
 (* True (at least for max=100) *)

I think this fulfils your "clean" and "non-contrived" criteria, but it is definitely not efficient: I set max to 100 because I didn't feel like waiting more than a few minutes for the answer!

EDIT

Also using Table, but without the If:

Clear[max, vals, a2];
max = 100;
vals = Table[2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 
  2 n (x + y + z - 3)), {x, 1, max}, {y, 1, x}, {z, 1, y}, {n, 1, 
max}];
a2 = BinCounts[Flatten@vals, {1, max + 1, 1}]

EDIT for SparseArray

Here is an approach using SparseArray in place of Table to get vals in the above. It is somewhat more efficient than Table, but not as efficient as the For loop way:

Clear[max, val, vals, a3];
max = 100;
vals = SparseArray[{x_, y_, z_, n_} /; 2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 
      2 n (x + y + z - 3)) <= max && z <= y <= x :> 2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 2 n (x + y + z - 3)), {max, max, max, max}];
a3 = Normal@BinCounts[Flatten@vals, {1, max + 1, 1}]

a3==a
 (* True *)

There is probably a way to make the condition more readable, but I haven't found it.

Consider the relative timings for max==100:

• For loops ~ 0.006 s
• Table ~ 162 s
• SparseArray ~ 0.8 s

But even SparseArray becomes horribly slow for n = 200.

Answer 2

Here's a version where the iterations are expressed recursively instead of imperatively:

Module[{val, iter, max, a}
, max = 5000
; a = ConstantArray[0, max]
; val[x_, y_, z_, n_] :=
    2 (2 n^2+(y-2) (z-2)+x (y+z-2)+2 n (x+y+z-3))
; iter[x_] /; val[x, 1, 1, 1] <= max :=
    (iter[x, 1]; iter[x + 1])
; iter[x_, y_] /; y <= x && val[x, y, 1, 1] <= max :=
    (iter[x, y, 1]; iter[x, y + 1])
; iter[x_, y_, z_] /; z <= y && val[x, y, z, 1] <= max :=
    (iter[x, y, z, 1]; iter[x, y, z + 1])
; iter[x_, y_, z_, n_] :=
    val[x, y, z, n] /. v_ /; v <= max :> (++a[[v]]; iter[x, y, z, n+1])
; Block[{$RecursionLimit = Infinity}, iter[1]]
; ListPlot[a]
] // Timing

On my machine, it runs about twice as slow as the For version.

Answer 3

I have a solution that I think is somewhat cleaner, and which still completes in a reasonable (but considerably longer) amount of time, and based on my desultory testing, it seems to scale with max at the same rate that the original version does. However, where on my machine the original version takes about 3 sec. to complete for max = 5000, my version takes about 40 sec. to complete.

Here's my solution.

valuesToCounts[vals_, max_] := Normal@SparseArray[Rule @@@ Tally[vals], {max}]

countN[fun_, max_, range_, x_, y_, z_] :=
 With[{closure = fun[x, y, z, #] &},
  valuesToCounts[closure /@ TakeWhile[range, closure@# <= max &], 
   max]]

countZ[fun_, max_, range_, x_, y_] :=
 Total[countN[fun, max, range, x, y, #] & /@ 
   TakeWhile[range, fun[x, y, #, 1] <= max && # <= y &]]

countY[fun_, max_, range_, x_] :=
 Total[countZ[fun, max, range, x, #] & /@ 
   TakeWhile[range, fun[x, #, 1, 1] <= max && # <= x &]]

countX[fun_, max_, range_] :=
 Total[countY[fun, max, range, #] & /@ 
  TakeWhile[range, fun[#, 1, 1, 1] <= max &]]

At the risk of belaboring the point, and compromising the functional purity of my solution, I tried to improve performance by using a closure to emulate pass-by-reference (a Mathematica trick I heartily recommend), like so:

scanN[fun_, scanner_, max_, range_, x_, y_, z_] :=
 With[{closure = fun[x, y, z, #] &},
  Scan[scanner, closure /@ TakeWhile[range, closure@# <= max &], max]];

scanZ[fun_, scanner_, max_, range_, x_, y_] :=
 scanN[fun, scanner, max, range, x, y, #] & /@ 
  TakeWhile[range, fun[x, y, #, 1] <= max && # <= y &];

scanY[fun_, scanner_, max_, range_, x_] :=
 scanZ[fun, scanner, max, range, x, #] & /@ 
  TakeWhile[range, fun[x, #, 1, 1] <= max && # <= x &];

scanX[fun_, scanner_, max_, range_] :=
 scanY[fun, scanner, max, range, #] & /@ 
  TakeWhile[range, fun[#, 1, 1, 1] <= max &];

countByScanning[fun_, max_] :=
  Module[{a = ConstantArray[0, max], range = Range[max]},
   scanX[fun, (a[[#]]++) &, max, range];
   a];

There's a fair amount of code repetition here; I suspect with a little more work it could be massaged into something even prettier. Still, it's a good deal slower, so it might not be worth the candle.

Answer 4

My solution is not elegant but at least it is rather fast. The idea is similar to the other answers. I create one big and clumsy iterator instead of four simple ones. In order to compare timings I need to say that AbsoluteTiming of the original code on my machine for max=5000 is 4.7806686.

ClearAll[next, step, val];
max = 5000;

pat = {x_, y_, z_, n_};
val[pat] := 
  2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 2 n (x + y + z - 3));
a1 = ConstantArray[0, max];

next[pat] := Which[
   x != y == z == n == 1, 0,
   (x >= y && z == n == 1) || x == y == z, {x + 1, 1, 1, 1},
   (x > y > z && n == 1) || x > y == z, {x, y + 1, 1, 1},
   x >= y > z, {x, y, z + 1, 1}
   ];

step[p : pat] := If[
   (r = val@p) <= max, a1[[r]]++; {x, y, z, n + 1},
   next@p
   ];

NestWhile[step, {1, 1, 1, 1}, ! (# === 0) &]; // AbsoluteTiming
a1 == a

(*==>   {7.2491184, Null} 
        True *)

The slow-down ratio is 1.5. The difference in speed here is because of ugly exit condition and redundant comparisons. As long as we want only to avoid For we can do something like this:

ClearAll[next, val];
ClearSystemCache[];
max = 5000;
x = y = z = n = 1;
val := 2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 
     2 n (x + y + z - 3));
a1 = ConstantArray[0, max];

run = True;

next := Which[
   x != y == z == n == 1, run = False,
   x == y == z || (x >= y && z == n == 1), (x++; y = z = n = 1),
   x > y == z || (x > y > z && n == 1), (y++; z = n = 1),
   x >= y > z, (z++; n = 1)
   ];

While[run,
  If[(r = val) <= max, a1[[r]]++; n++, next]
  ] // AbsoluteTiming
a1 == a

(* {5.1868692, Null}
   True *)

Nevertheless I did not manage to beat the original code in speed (I don't speak about the elegance for obvious reasons).