# Can this be written well, without loops?

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.

-
Oh dear........ –  rm -rf Feb 22 '12 at 19:52
Can you give any more info as to where this code came from and how it is intended to be used? –  rcollyer Feb 22 '12 at 20:01
@rcollyer I'm sorry, I cannot. Functionally speaking it's all right there and fairly straightforward once you get past the ugliness. –  Mr.Wizard Feb 22 '12 at 20:04
@rcollyer Most likely a ProjectEuler problem. For some strange reason, ListPlot[a] looks very familiar... –  rm -rf Feb 22 '12 at 20:06
I have a very fast solution that only works for max=5000. It involves DumpSave and Get –  Rojo Feb 22 '12 at 21:56

There are some features of this specific problem one can take advantage of. The boundary of the x,y,z,n domain represented by val <= max is linear in x,y,z and only quadratic in n; furthermore val increases with each of the variables. So basically the loops might be done in any order, and the limits might be solved for explicitly.

We'll start with the limit max and the expression val, which can be compiled for the sake of comparison.

max = 5000;
val[x_, y_, z_, n_] :=
2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 2 n (x + y + z - 3));
valc = Compile[{{x, _Integer}, {y, _Integer}, {z, _Integer}, {n, _Integer}},
2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 2 n (x + y + z - 3))];


We can then solve for the limits on the indices z,y,x,n and save them in idxLimit[tag], where tag runs 1 through 4 and corresponds to z,y,x,n in that order. (Here a function runs through the tags and sets up idxLimit, but it could have been set up with separate formulas just as easily, as in the output below the code.)

Function[{tag}, idxLimit[tag] =
Simplify[
Min[{y, x, {}, {}}[[tag]], {z, y, x, n}[[tag]] /.
Last@Solve[(val[x, y, z, n] /. Take[{z -> 1, y -> 1, x -> 1}, tag - 1]) == max,
{z, y, x, n}[[tag]] ]],
n >= 1]
] /@ Range[4];


Next we make the table of the values

a = Normal@SparseArray[Rule @@@ #] &@
Tally[Flatten[
With[{i1 = {z, idxLimit[1]},
i2 = {y, idxLimit[2]},
i3 = {x, idxLimit[3]},
i4 = {n, idxLimit[4]}},
Table[valc @@ {x, y, z, n}, i4, i3, i2, i1] ] ]
]; // AbsoluteTiming
(* {1.536007, Null} *)


Below is a table of timings (in sec.) that compares using val instead of valc and ParallelTable (on a 2-core machine). It also compares the timing of the OP's For-loop program, with and without a compiled val. The last line are the timings for max = 20000.

Table itself accounts for about 0.876687 sec. (10.003417 for 20K). Most of the rest of the time is for evaluation valc (about 7-8 sec. in the 20K case) or val. A smaller chunk is spent collecting the results. It seemed while I was playing with the problem, that a[[r]]++ suffers from having to evaluate a[[r]] twice, once for reading and once for writing. Perhaps it doesn't take that much time, but I felt like there was a limit to how fast I could accumulate results in a that way. The Table way gains a little time at the expense of a lot more memory.

It's debatable whether solving for the limits (to get idxLimit) is clean. The original val leads to the strange expressions. Mainly it's a mathematical trick than a programming one, which allows a rather standard conversion of for loops to Table. The rest of it are just tweaks.

-
Thank you, Michael. –  Mr.Wizard Jan 22 '13 at 11:08
I am going to Accept this. Though it is not a purely programmatic solution and it didn't teach me a new coding trick as Simon's answer did it does finally provide a faster method than my For loops, and leverages Mathematica's strengths to do so. –  Mr.Wizard Feb 9 '13 at 11:51

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.

-
Yes, the problem here is that we do not exit the loop when val > max -- one can hack it with Return or the like but it makes it more ugly instead of less. –  Mr.Wizard Feb 22 '12 at 21:43
The problem (apart from efficiency issues you noted) is that, in general, you can not know that x, y, z and n should only change within the interval {1,max}. You can probably prove it in some cases, but that's probably besides the point. –  Leonid Shifrin Feb 22 '12 at 21:44
You can always use Break in your table to escape early. I wouldn't necessarily call that clean, though! :) –  Pillsy Feb 22 '12 at 21:47
@Pillsy right, but then you are better off using For. –  Mr.Wizard Feb 22 '12 at 21:48
@Leonid I do hope you have a solution for me. –  Mr.Wizard Feb 22 '12 at 21:48

I cannot write the code without any loops, but I have pared it down to a single While loop, making extensive use of short-circuiting to bail out of the expression at the appropriate point. Whether you consider this cleaner than nested For loops is a matter of taste. I personally find it easier to follow.

max=5000;
a={};
x=y=z=n=1;
val:=2 (2 n^2+(y-2) (z-2)+x (y+z-2)+2 n (x+y+z-3));

While@Or[
((r=val)<=max)&&(n++;a={a,r};True),
(n=1;++z<=y&&val<=max),
(z=1;++y<=x&&val<=max),
(y=1;x++;val<=max)
];

a=Tally@Flatten@a;


This gives a as a list of {position,value} pairs which you could use with something like ListPlot, or convert to the original form with a = Normal@SparseArray[Rule@@@a]

The code is a smidge faster than the original, but this is entirely due to building a as a linked list. The nested For loops are even faster if you make the same change.

For whatever reason the code above crashes Mathematica 7 when max is higher.
Here is working v7 code using pre-allocation:

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));

While @ Or[
(r = val) <= max && (n++; a[[r]]++; True),
n = 1; ++z <= y && val <= max,
z = 1; ++y <= x && val <= max,
y = 1; x++; val <= max
];

-
That looks impressive! I need time to digest this. –  Mr.Wizard Jul 3 '12 at 2:15
I must second that impression. That is some seriously interesting code. –  rcollyer Jul 3 '12 at 2:27
I'm having trouble separating the linked list aspect of this from the underlying algorithm. I quickly observe that the memory requirements of this are much greater than my original, and more troubling Mathematica crashes if max is greater than about 7500. By the way ((r=val)<=max)&&(n++;a={a,r};True) is a very interesting fragment of code; I don't think I've ever used And for conditional evaluation. –  Mr.Wizard Jul 3 '12 at 8:53
@Mr.Wizard, the increased memory usage is because of the linked list approach. In your code a is a Packed Array, in mine it isn't - I simply accumulate values of r and tally them afterwards. You can replace a={a,r} in my While loop by a[[r]]++ to make it build a the same way your code does. I'm not getting any crashes. With max=10000 I get a ByteCount for a of about 80MB (compared to 40kB for your code) and MemoryInUse is about 100MB. –  Simon Woods Jul 3 '12 at 20:27
Re the compiling problem, it seems like there might be a difference between v7 and v8. This works on v8: f=Module[{x,y,z,n,a,r},With[{val=2(2n^2+(y-2)(z-2)+x(y+z-2)+2n(x+y+z-3))},Compi‌​le[{{max,_Integer}},Block[{a=Table[0,{max}],x=1,y=1,z=1,n=1,r=0},While@Or[(r=val)‌​<=max&&(n++;a[[r]]++;True),n=1;++z<=y&&val<=max,z=1;++y<=x&&val<=max,y=1;x++;val<‌​=max];a]]]]; a=f[5000]; –  Simon Woods Jul 4 '12 at 10:11

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.

-
This is nice: it makes the tree structure of the problem manifest. +1 –  JxB Feb 24 '12 at 0:22

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.

-

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).

-