How to implement this simple procedural style algorithm efficiently in Mathematica style code?

Question

I am trying to move my C solution code for project euler problem 135 to Mathematica. But I encountered serious performance problem. The algorithm is very simple, here is a plain translation:

sol[_] := 0
For[x = 1, x < 750000, x++,
  For[d = Ceiling[x/3], d < 250000, d++,
   n = 3*d*d + 2*d*x - x*x;
   If[n < 1000000, sol[n]++, d = 250000]
   ]
  ];
Count[sol /@ Range[1000000], 10]

Of course this is not functional style and should be considered low efficiency since it uses procedural loops (this code takes 13s on my laptop). However, it turns out that the more I rewrite this algorithm in functional style, the slower it is.

One other thing worth mentioning: I tried to compile the code by changing the downsvalues to a list of length 1000000. When I checked with CompilePrint, everything was fine and there is no MainEvaluate, but messages saying: non-tensor object is generated in runtime, and I have no idea why it is generated.

My question is: is it possible to implement the same/similar algorithm in good Mathematica style, and have reasonable performance?

Appendix

Here is the Library Link code for the same algorithm, runs in 0.006s on my laptop. I don't expect the Mathematica code to be as fast as this, but as least within several seconds...

src = "#include \"WolframLibrary.h\"
  #define X 1000000
  DLLEXPORT int euler135(WolframLibraryData libData, mint Argc, \
MArgument *Args, MArgument Res) {
      mint num_of_sol, x_limit, d_limit;
      mint result = 0;
      x_limit = MArgument_getInteger(Args[0]);
      d_limit = MArgument_getInteger(Args[1]);
      num_of_sol = MArgument_getInteger(Args[2]);
      
      int sol[X]={0};
      mint n=0;
      for (mint x=1; x<x_limit; x++) {
          for (mint d=x/3+1; d<d_limit; d++) {
              n=3*d*d+2*d*x-x*x;
              if (n<X) {
                  sol[n]++;
              }else{
                  d=d_limit;
              }
          }
      }
      for (mint i = 1; i<X; i++) {
          if (sol[i]==num_of_sol) {
              result++;
          }
      }
      MArgument_setInteger(Res,result);
      return LIBRARY_NO_ERROR;
  }";
Needs["CCompilerDriver`"];
lib = CreateLibrary[src, "euler135"];
euler135 = 
  LibraryFunctionLoad[lib, "euler135", {Integer, Integer, Integer}, 
   Integer];
euler135[750000, 250000, 10] // AbsoluteTiming

Answer 1

You probably have something like this by now but just to show, something of this sort can be put through Compile. However: (1) I believe this is not quite what is needed to hadnle PE 135 (2) It is not counting the time needed to create the CompiledFunction, hence might be regarded as cheating. (3) I'm probably not giving the right variable to compile, my guess is it should be based on the one million mentiond in PE 135 and not on 750000 which I gather is derived by analysis.

countSol = Compile[{{xmax, _Integer}, {nsols, _Integer}}, Module[
    {sol = ConstantArray[0, 10^6 - 1], n, dmax = Ceiling[xmax/3]},
    Do[
     Do[n = 3*d^2 + 2*d*x - x^2;
      If[0 < n < 10^6, sol[[n]]++;, Break[]],
      {d, Ceiling[x/3], dmax}], {x, xmax}];
    Count[sol, nsols]], CompilationTarget -> "C", 
   RuntimeOptions -> "Speed"];

AbsoluteTiming[
 countSol[750000, 10]
 ]

(* {0.034309, 4290} *)

Answer 2

Your translation runs in 40.3125 on my computer.

The following modification runs in 10.7031.

 Count[Tally[Join @@ MapThread[Table[d (3 d + 2 #) - #2, {d, #3, #4}] &,
           {#, #^2, Ceiling[#/3], Floor[(2 Sqrt[750000 + #^2] - #)/3] &[N[#]]} &[
          Range[750000 - 1]]]][[All, 2]], 10]

You square x for each d, which is not needed, so I have put it outside the inner loop. Also you use many top-level Ceiling-calls, which I have replaced by a single Ceiling with a list argument. I don't think performance is significantly affected by the fact that I use arithmetic to get the upper d value rather than testing in the inner loop.

As to reasonable performance I think Compile is necessary.