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
n>0
. $\endgroup$