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++,
   ]; n = 1
  ]; z = 1
 ]; y = 1

The output is the array a.

  • 20
    $\begingroup$ Oh dear... ${}$ $\endgroup$
    – rm -rf
    Commented Feb 22, 2012 at 19:52
  • 2
    $\begingroup$ Can you give any more info as to where this code came from and how it is intended to be used? $\endgroup$
    – rcollyer
    Commented Feb 22, 2012 at 20:01
  • 1
    $\begingroup$ @David I cannot go into a lot of detail about the problem. Consider val a black-box function and just look at the loop structure. I think what that does it reasonably clear. $\endgroup$
    – Mr.Wizard
    Commented Feb 22, 2012 at 20:40
  • 7
    $\begingroup$ I have a very fast solution that only works for max=5000. It involves DumpSave and Get $\endgroup$
    – Rojo
    Commented Feb 22, 2012 at 21:56
  • 4
    $\begingroup$ @Rojo, your comment reminds me of a beautiful remark by O'Keefe in The craft of Prolog. It's in a discussion of the problem: find a 9-digit number, with distinct digits, such that the first n digits form a number divisible by n, for 1 ≤ n ≤ 9. He goes through several optimisations, but concludes with (paraphrased) "you have to be careful, or eventually you will wind up optimising your program to :- digits(381654729)."(en.wikipedia.org/wiki/Polydivisible_number#Background). $\endgroup$
    – LSpice
    Commented Apr 23, 2015 at 23:47

8 Answers 8


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] =
     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];

info on idxLimit

Next we make the table of the values

a = Normal@SparseArray[Rule @@@ #] &@
      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 of timings

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.

  • $\begingroup$ 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. $\endgroup$
    – Mr.Wizard
    Commented Feb 9, 2013 at 11:51

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.

val:=2 (2 n^2+(y-2) (z-2)+x (y+z-2)+2 n (x+y+z-3));



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
  • 1
    $\begingroup$ That looks impressive! I need time to digest this. $\endgroup$
    – Mr.Wizard
    Commented Jul 3, 2012 at 2:15
  • $\begingroup$ I must second that impression. That is some seriously interesting code. $\endgroup$
    – rcollyer
    Commented Jul 3, 2012 at 2:27
  • $\begingroup$ 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. $\endgroup$
    – Mr.Wizard
    Commented Jul 3, 2012 at 8:53
  • $\begingroup$ @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. $\endgroup$ Commented Jul 3, 2012 at 20:27
  • 1
    $\begingroup$ 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))},Compile[{{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]; $\endgroup$ Commented Jul 4, 2012 at 10:11

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}];

 (* 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!


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, 
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}]

 (* 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.

  • $\begingroup$ 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. $\endgroup$
    – Mr.Wizard
    Commented Feb 22, 2012 at 21:43
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Feb 22, 2012 at 21:44
  • $\begingroup$ You can always use Break in your table to escape early. I wouldn't necessarily call that clean, though! :) $\endgroup$
    – Pillsy
    Commented Feb 22, 2012 at 21:47
  • $\begingroup$ @Pillsy right, but then you are better off using For. $\endgroup$
    – Mr.Wizard
    Commented Feb 22, 2012 at 21:48
  • $\begingroup$ @Leonid I do hope you have a solution for me. $\endgroup$
    – Mr.Wizard
    Commented Feb 22, 2012 at 21:48

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.

  • $\begingroup$ This is nice: it makes the tree structure of the problem manifest. +1 $\endgroup$
    – JxB
    Commented Feb 24, 2012 at 0:22

This is an interesting, if (necessarily) ill-posed question. My approach is to couch it in more general terms and attempt to clarify what is possible, natural and generalizable by exploring the meaning of "clean", "efficient" and "non-contrived". Some implications for language design are also discussed.

max = 5000;
val := 2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 2 n (x + y + z - 3));
ai[i_] := Length@Solve[val == i && x >= y >= z >= 1 && n >= 1, {x, y, z, n}, Integers];

a = Table[ai[i], {i, 1, max}]

   - Warning: On a personal machine this is likely to take several hours;
   for plausibility setting max = 100 takes less than 1 min to complete.

   - Acknowledgment: This was computed at the Pawsey Supercomputing Center
    ( ~18 min on 12 kernels with ParallelTable replacing Table).


Hence as per the OP's request, the code is written without loops and is also, I'd argue, clean, non-contrived, (space) efficient; offers lazy evaluation for a's components and provides additional context. It is however, manifestly and grossly inferior from a time-efficiency viewpoint needing HPC resources to even confirm its equivalence with previous implementations. Its conceptualization however, dilineates likely limitations.


One point alluded to, but not emphasised in the answers/comments is that constructs like Table, Array, SparseArray etc. all "explicitly define" a space, whereas For, While etc. loops describe processes via stopping conditions. The former commonly:

  • Specify a space of elements
  • Specify how that space is traversed
  • Specify computations on visited elements during the traversal

directly and more naturally (arguably because humans are more visual than logical creatures). Looping constructs like For, on the other hand, can be co-opted into defining spaces with the above features but often do so sub-optimally as their stopping conditions end up getting in the way.

The underlying space being defined in the question is not immediately obvious, as evidenced by initially resorting to loops but the deeper question seemingly being proposed is how permanent is this state of affairs? is there a systematic way of re-casting into more functional-like Table/Array forms? Unfortunately the answer turns out to be no, although there is usually more that can be done; first though, consider how this particular loop could be translated into a Table form albeit in a way apparently insufficiently "magical" or "programmatical" - by analyzing val's particular structure and mathematical properties.

First, it can be noted that the For loop terminates since val's quadratic increase in n is guaranteed to eventually exceed max. The actual iterates over which this takes place however, is not immediately obvious corresponding to the non-obviousness of converting into a Table/Array formulation. MichaelS2's answer, alone amongst all the responses, explicitly finds iterates thereby arriving at a non-loop solution as originally requested (it is on this basis, I'd argue, that his answer merited acceptance rather than from the stated efficiency gains (other Table/SparseArray solutions assumed, stopping properties related to val [as acknowledged in the comments] while all the other answers have loops lurking somewhere within While's or recursive procedures).

[Note that even this loop dichotomy can't be pushed too far; a stopping condition can be inserted into a Table's iterates thereby turning it into a "Loop" whilst explicit iterates in a For loop can augment stopping-conditions thereby turning it into a "Table"]

Converting to a Table formulation by finding explicit iterates however, required analyzing val in the stopping condition with the specific nature of such analysis unavoidable thereby dashing the OP's hopes and intuition (admittedly also my own) for a generic, "magical", "programming solution" for loop conversion. This follows since to assume otherwise would mean being able to translate arbitrary For stopping conditions into decidable procedures (those characteristic of Table-like iterates), something impossible from the Halting Problem's unsolvability. For some examples therefore, the traversal simply has to be stepped through (imagine a random val to simulate a black box function).

The Solve reformulation conceptualizes a as a vector who's i'th component counts the number of "essentially different" solutions to the Diophantine equation: val == i. Here "essentially different", means up to re-labeling the x, y,z variables and is implemented by translating For's stopping conditions that effectively insist on the variable ordering x>=y>=z>=1.

The explict values of Table iterates correspond to upper bounds on the size of integers in any putative solution. Diophantine equations are often used to settle decidability questions which here emerges in the form of systematising code re-writing attempts.

In terms of (time) efficiency implications, one pertinent issue is whether or not these counts can be performed without explicitly generating actual solutions. SatisfiabilityCount offers an interface for doing so in relation to Boolean equations although its performance suggests solutions are first generated before being counted.

A Demonstration shows that while it is possible to generate closed-form formulae for counting solutions without their explicit generation for classes of equations (albeit over small dimensions).

In general however, counting problems of #P complexity are generally intractable and the point of couching the problem in number-theoretic terms is that it can show you what you are up against in terms of searching for efficiency gains. I suspect that existing complexity results for counting solutions of diophantine equations makes the prospects of significant improvements in this example very limited.

Intuitively, a depends on visiting every element in the For loop with any significant efficiency gains stemming from being able to short-cut this process. The exhibited efficiency gain (with the Table solution) does this in a limited (if clever and useful) way by using max (effectively bypassing the stopping condition as the means for ignoring those variable values for which val exceeds max) while also exploiting Mathematica's implementation of Table (its compilability and parallelization).

The core "irreducibility" of this computation however, can be first discerned by observing a's "randomness" for the first 5K elements:

enter image description here

with some order emerging from seeing the first 50K elements:

enter image description here

but all the while random-like upper boundary persisting as evident from viewing the first 250K elements:

enter image description here

(* Acknowledgement:
   - The above plots used the compiled "Table" from Michael's answer.
   - The last 250K plot used the resources of the Pawsey Supercomputing Centre
     in Perth, Western Australia (taking ~15 min with 12 Kernels and ParallelTable) 

indicating the absence of a recursive reducibility (e.g. impossibility of expressing a[[i]] in terms of a[[j]]'s for j<i ). Consider however, a similar example instead involving the expression val2 (derived from val by replacing n^2 with n and removing the later occurrence of 2n). Now there is evidence of clear reducibility.

enter image description here

The reducibility inherent in the the "val2 computation" suggests opportunities for efficiency gains in codified mathematical knowledge, say that built-in to functions like Solve; sure enough it outperforms the For loop for finding a[[5000]] (17.73 s vs 120.27 s) in contrast to its inferior performance in relation to the more "irreducible val computation" (5.82 s vs 2.07 s).

enter image description here

To be sure, finding a single element of a is Solve's focus in contrast to For's focus here in generating all of a's elements. On the other hand, this also indicates one might have anticipated Solve's more competitive performance (notwithstanding its impressive generality and suggesting its improvement by associating predicates corresponding to For's stopping conditions together with preliminary checks on identifying/looking-up irreducibility); at any rate, the same effects would be observable given a mature counting framework in Mathematica (i.e. functions for which counting was the focus).

There are also many implications for language design here but this post is plenty wordy already.

Well some implications:

The conventional wisdom seems to be that For's should be afforded sideways glances before being cast off into outer darkness while moving into the utopia of Table-Array vectorization. While this view carries a certain force (I can't imagine programming without Tables/Arrays) IMO such unyielding focus on functional programming can also become limiting.

Firstly, while Tables/Arrays represent powerful ways of exploring the computational universe, they do so in a regimented way progressively fixing variable dimensions that while easily humanly graspable potentially ignore potentially fruitful search spaces. In addition, they also tend to push answers in directions that may not even require brute-force enumeration. Even in situations where brute-force enumeration is apparently unavoidable however, (such as the irreducibility illustrated here) certain types of questions necessitate a different, "more semantic" approach. This is possibly foreshadowed in pure math initiatives but more broadly it also has implications in the scientific practice of model-building.

The process of identifying this problem (from apparently Project Euler) involved a type of reverse-engineering, analogs of which frequently occur in general modeling. To relate it back to the problem at hand: It's not necessarily the case that saying (or computing) something useful about a[[i]] should require computing its exact value. As a toy example, in a more "semantic" computing environment, ZeroQ[a[[(googol=10^100)-1]] should return True without setting out to generate the structure a or the exact value of a[[googol-1]] (n.b. current output of PossibleZeroQ[googol] or PossibleZeroQ[googol^googol])

This suggests notions of "LazyTable", "LazyArray" analogs (SemanticArray, DelayedArray or Array overloaded?) that move gracefully between producing actual structures when immediately, computationally feasible and otherwise; accessing pre-computed databases and/or remaining unevaluated in readiness for input into computational questions about a[[i]] (i.e. that don't require its exact value).

While the ZeroQ question/computation here is toy, situations inevitably arise in which the computation forms a link in a useful (computational/reasoning) chain and/or it is in fact the only way to identify the origin of the original computation. Consider a (hidden) For loop representing an irreducible natural process in which say a[[googol;;googol+100]], constitutes the output of 100 experimental observations (i.e. a's index corresponds to time). Due to irreducibility it may be impossible to run the computation again to deduce such provenance but nonetheless, computations about a[googol;;googol+100] (i.e. beyond ZeroQ) potentially exhibit identifying fingerprints traceable back to the originating For process. This however, requires For's original framing and collation (or equivalents) in "lazy" terms.

This computational gap is a fundamental limitations of the (still) amazingly useful look-up tables such as the On-line Encyclopadia of Integer Sequences, (n.b. FindSequenceFunction) but functions like "LazyTable" at least provide a mechanism for overcoming such gaps by generalizing this collation. It portends the usefulness of (inter-related) look-up tables for lists of functions and/or symbols backed by their networked connections (integers as algorithmic fingerprints owe as much to a human's proclivity for counting things as anything else); in so doing significant parts of scientific inquiry can be automated given how modeling/simulations often embody precisely this reverse-engineering process.

  • $\begingroup$ Thanks for this answer. It's going to take some time to fully appreciate but I already find it of value. And you note that I did Accept Michael's answer and not simply on the grounds of speed but on the analytical approach taken. $\endgroup$
    – Mr.Wizard
    Commented May 18, 2015 at 16:53
  • $\begingroup$ I find your closing line tantalizing: "There are also many implications for language design here but this post is plenty wordy already." $\endgroup$
    – Mr.Wizard
    Commented May 18, 2015 at 16:55
  • $\begingroup$ Yes, and of course irrespective, it is the OP's perogative :) $\endgroup$ Commented May 18, 2015 at 16:55
  • $\begingroup$ By the way, and no offense intended, this kind of direct Solve formulation was what I was trying to avoid with the requirement of efficiency. I would very much have liked it to be fast but it was not and I eventually aborted it. At the time I posted this I was very much against posting Project Euler problems, hence the obfuscation, but I've since come to realize it is a losing battle, so I suppose the question might be clarified. $\endgroup$
    – Mr.Wizard
    Commented May 18, 2015 at 16:59
  • $\begingroup$ None taken and yes I did suspect such reverse-engineering might be a possibility and I would agree that it can't hurt to be up-front. Doesn't this however, say something about Solve's implementation - shouldn't its efficiency be comparable to a (For) "brute force" enumeration? $\endgroup$ Commented May 18, 2015 at 17:13

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 &], 

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

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},

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

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


a 'loopless' one liner.. This takes ~5 minutes, way slower than the original, but considerably faster that Ronald's..

 max = 5000;
 Clear[val, x, y, z, n]
 val[x_, y_, z_, n_] := 
      2 (2 n^2 + (y - 2) (z - 2) + x (y + z - 2) + 2 n (x + y + z - 3));
 a = Normal@
    SparseArray[Rule @@@ #, max] &@(val[x, y, z, n] /. 
          FindInstance[  val[x, y, z, n] <= max &&  
              x >= y >= z >= 1 && n >= 1, {x, y, z, n} , Integers, 10^6] //

Unfortunately we need to supply FindInstance with a bound on the number of instances..

  • $\begingroup$ Indeed, but the bound means the computation is "not certain" - some certainty has been sacrificed for speed. Admittedly though, it doesn't seem like much certainty has been sacrificed looking at the graphs. The 1M bound seems overly conservative although curiously, dropping it seems to add a few minutes to the computation. $\endgroup$ Commented May 25, 2015 at 9:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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