How to speed up my Project Euler code

The evaluation speed of Mathematica often depresses me. It did it again when I wrote a code to solve a Project Euler problem. https://projecteuler.net/problem=206

Here's my code:

ClearAll[a];
a =
Compile[{},
NestWhile[# + 1 &, 10^9, ! MatchQ[
IntegerDigits[#^2], {1, _, 2, _, 3, _, 4, _, 5, _, 6, _, 7, _,8, _, 9, _, 0}] &],
CompilationTarget -> "C"]


I think it is the simplest solution. I can't find other solutions which are simpler. But my Mathematica runs a long time before getting the answer. How can I speed this up? I've run into similar problems while I was solving other problems.

• Pattern matching is often slow, you might try something like IntegerDigits[#^2][[ ;; ;; 2]] == {1,2,3,4,5,6,7,8,9,0} instead. – C. E. Sep 3 '14 at 7:16
• #+1 can be replaced by #+10, because a square number ends with "0" must ends with "00" – wuyingddg Sep 3 '14 at 8:20
• If you realize that the maximum solution is 10 Floor[Ceiling[Sqrt[1929394959697989990]]/10] and step down by 10 from there, you'll get the solution almost immediately. – Mark McClure Sep 3 '14 at 12:38
• @MarkMcClure Yes but this is because you know the answer. It could also be very close to 102030405060708090. – anderstood Oct 21 '15 at 20:00

Compilation is a certainly a good idea if you're going the brute-force route. So let's first tackle that, which will get us the answer in roughly 30 seconds computation time. Afterwards we'll come up with better strategy, which can do it in 0.3 seconds.

Compilation

A performance pitfall when compiling functions is that sometimes the compiled function just calls the main evaluation loop, effectively gaining nothing. We can check this with CompilePrint:

Needs["CompiledFunctionTools"]

a // CompilePrint

      (...)
Result = I3

1    I3 = MainEvaluate[(...)]
2    Return


The call to MainEvaluate is the culprit here.

So let's rewrite your function such that it does compile properly. Using Pickett's and wuyingddg suggestions, we end up with:

PerformSearch = Compile[
{
{startValue, _Integer},
{increment,  _Integer}
},
NestWhile[
# + increment &,
startValue,
IntegerDigits[#^2][[-1 ;; 1 ;; -2]] =!= {9, 8, 7, 6, 5, 4, 3, 2, 1} &
]
]


Let's check if this did compile ok:

PerformSearch // CompilePrint

      (...)

1    I4 = I0
2    I3 = I4
3    I7 = Square[ I3]
4    T(I1)2 = IntegerDigits[ I7, I5]]
5    T(I1)0 = Part[ T(I1)2, T(I1)1Span]
6    B0 = CompareTensor[ I6, R2, T(I1)0, T(I1)3]]
7    if[ !B0] goto 12
8    I3 = I4
9    I7 = I3 + I1
10   I4 = I7
11   goto 2
12   Return


Ok, that looks much better. We can now perform the search, but not after we've determined the range of possible solutions:

max =   Floor @ Sqrt @ FromDigits @ Riffle[Range[9], 9] (* 138902662 *)
min = Ceiling @ Sqrt @ FromDigits @ Riffle[Range[9], 0] (* 101010102 *)


As an aside, note that max^2 is below $MaxMachineInteger on 64-bit systems. But on 32-bit it isn't, which causes PerformSearch to switch back to the uncompiled code. Keeping Mark McClure's comment in mind, we'll cheat a bit and start from the maximum to find immediately: result = PerformSearch[max, -1]  138901917  How much faster is this than starting from the minimum? (result - min)/(max - result) // N  50861.5  Roughly 50 thousand times faster, not bad! Starting from the minimum takes about 30 seconds on my machine. Last but not least, let's double-check if we've obtained the correct number: result^2  19293742546274889  If you multiply this with 100, you've got your desired integer. Non-brute-force The approach above scanned roughly 38 million (!) integers in the worst-case scenario (starting from the minimum). Other answers to the OP's question have shown that you can and should go one better than that. Here's my take on an efficient general solution, using an iterative step-by-step process: ClearAll[ FindIntegerRoots, CheckIfEmpty, SelectDigitSequences, InsertNewDigits, ConvertPattern ]; FindIntegerRoots[ pattern : { (0|1|2|3|4|5|6|7|8|9|Verbatim[_]).. }, power_Integer: 2 ] /; power > 1 := With[ { maxRoot = Floor[FromDigits[pattern /. Verbatim[_] -> 9]^(1/power)] }, FromDigits /@ SelectDigitSequences[ Fold[ CheckIfEmpty @ SelectDigitSequences[ InsertNewDigits @ #1, maxRoot, power, ConvertPattern[pattern, #2] ] &, {{}}, Range @ IntegerLength @ maxRoot ], maxRoot, power, ConvertPattern[pattern, Length @ pattern] ] ~Quiet~ {CompiledFunction::cfnlts} ~Catch~ "isEmpty" ]; CheckIfEmpty[{}] := Throw[{}, "isEmpty"]; CheckIfEmpty[nonEmpty_] := nonEmpty; SelectDigitSequences = Compile[ { {digitSequences, _Integer, 2}, {maxRoot, _Integer }, {power, _Integer }, {numDigits, _Integer }, {positions, _Integer, 1}, {compareDigits, _Integer, 1} }, Module[ { powersOfTen = 10^# & /@ Range[Length @ First @ digitSequences - 1, 0, -1] }, Select[ digitSequences, And[ # <= maxRoot, IntegerDigits[#^power, 10, numDigits][[positions]] == compareDigits ] & @ ( powersOfTen.# ) & ] ] ]; InsertNewDigits = With[ { range = List /@ Range[0, 9] }, Compile[ { {digitSequences, _Integer, 2} }, Outer[#1 ~Join~ #2 &, range, digitSequences, 1] ~Flatten~ 1 ] ]; ConvertPattern[pattern_, length_] := With[ { trimmed = pattern[[Length@pattern - length + 1 ;; -1]] }, Sequence @@ { length, Flatten @ Position[trimmed, _Integer], DeleteCases[trimmed, Verbatim[_]] } ];  While the code might be long, it is nevertheless fast: FindIntegerRoots @ Riffle[Range[9], _] // AbsoluteTiming  {0.284622, {138901917}}  A performance increase of 100 over the above brute-force method! But how efficient is it? Here's a nice graph of the amount of checks it performs in this case: That's about 296 thousand checks in total; certainly much less than 38 million! Some more concealed squares The nice thing about FindIntegerRoots is that it works on any pattern, not just the OP's case: FindIntegerRoots @ {_, 9}  {3, 7}  FindIntegerRoots @ { 1, _, 4}  {12}  FindIntegerRoots @ { 1, _, 6}  {14}  And if you're adventurous, you may even ask for roots of different powers: FindIntegerRoots[{_, _, _, 5}, 3]  {5, 15}  So let's have a bit of fun and search for more concealed squares. First, are there any partial concealed squares (i.e. of the form 1_2, 1_2_3, 1_2_3_4, etc)? FindIntegerRoots @ Riffle[Range @ #, _] & /@ Range @ 9  {{1}, {}, {}, {1312}, {}, {115256, 127334, 135254}, {}, {}, {138901917}}  Flatten[%]^2   {1, 1721344, 13283945536, 16213947556, 18293644516, 19293742546274889}  How about reverse concealed squares (i.e. 2_1, 3_2_1, 4_3_2_1, etc)? FindIntegerRoots @ Riffle[Reverse @ Range @ #, _] & /@ Range @ 9  {{1}, {}, {}, {}, {24171}, {}, {2657351, 2713399}, {}, {}}  Flatten[%]^2  {1, 584237241, 7061514337201, 7362534133201}  And finally, are there squares of the form 1_1_(...)_1_1? FindIntegerRoots @ Riffle[ConstantArray[1, #], _] & /@ Range[2, 9]  {{11}, {119, 131}, {}, {}, {110369}, {}, {10065739}, {}}  Flatten[%]^2  {121, 14161, 17161, 12181316161, 101319101616121}  It seems there are many nice squares. If you happen to find a particularly nice one, let me know :) • Which version of Mathematica are you using? It's quite strange that your code generates the CompiledFunction::cfn warning in my v8 (Win 32bit & 64bit) and v9 (Win 32bit), but it does work without warning on Wolfram Cloud! – xzczd Sep 3 '14 at 14:15 • @xzczd I've checked this on v9 on a Mac. The issue here is 32 vs 64 bit; $MaxMachineInteger should namely be bigger than FromDigits@Riffle[Range[9], 9], which it isn't on 32 bit. – Teake Nutma Sep 3 '14 at 16:06
• That seems to be the reason, BTW the \$MaxMachineInteger in my v8 (Win 64 bit) is the same as that in 32 bit, not sure if it's just the nature of v8… – xzczd Sep 4 '14 at 1:56
• Is using Big O notation in "Roughly O(10^4) times faster" a correct statement? It is surely 10^4 times faster in this example, but what exactly does O(10^4) times faster mean? – user Sep 4 '14 at 15:29
• @bruce14 What I meant here is that it's roughly 4 orders of magnitude faster; faster than 10^3 and slower than 10^5. – Teake Nutma Sep 4 '14 at 17:15

This non-paralell strategy takes less than 4 seconds in my stone age laptop (and doesn't use the "reverse searching" trick)-

First we check for the possible 6 digits endings (after discounting the final "00" ending in the squared number):

set= 2 Position[IntegerDigits[Range[10^5 + 1, 10^6-1, 2]^2][[All, -5;; -1;; 2]], {7,8,9}]+ 10^5 - 1;


Then we complete the whole set of numbers to check by adding all the possible "heads" between the max and min:

fullSet = Total[Tuples[{Range[101, 138] 10^6, Flatten@set}], {2}];


And then a quick search gets our number:

Cases[IntegerDigits[fullSet^2], {1, _, 2, _, 3, _, 4, _, 5, _,  6, _, 7, _, 8, _, 9}]

(*{{1, 9, 2, 9, 3, 7, 4, 2, 5, 4, 6, 2, 7, 4, 8, 8, 9}}*)


(* 1929374254627488900 *)


Edit

We can still cut the time in half if we see that the number should end in 3 or 7 for its square to end in 9:

set = Join[10 (Position[IntegerDigits[Range[10^5 + #, 10^6, 10]^2][[All, -5 ;; -1 ;; 2]],
{7, 8, 9}] - 1) + # & /@ {3, 7}] + 10^5;
fullSet = Total[Tuples[{Range[101, 138] 10^6, Flatten@set}], {2}];
Cases[IntegerDigits[fullSet^2], {1, _, 2, _, 3, _, 4, _, 5, _,  6, _, 7, _, 8, _, 9}]

• +1 Very concise and fast (~1.3 times faster than mine!). – Teake Nutma Sep 4 '14 at 20:31
• @TeakeNutma Thanks! I halved the time in the edit.Less than 2 secs now on my poor man's machine :) – Dr. belisarius Sep 5 '14 at 2:41
• That's pretty fast. But alas, I've updated my answer and now it's 2x faster than yours :). – Teake Nutma Sep 5 '14 at 10:11

There is a method by using Catch and Throw, and it gives result in about 90s. I think it can be improve a lot by Paralilize or ParallelEvaluate, but failed to finish that.

i = 10^9 + {30, 70};
isMatchQ =
If[IntegerDigits[#^2][[1 ;; -1 ;; 2]] ==  {1, 2, 3, 4, 5, 6, 7, 8,
9, 0}, Throw[#]] &;
Catch[Do[isMatchQ /@ i; i += 100, {10^9}]]


Edit

I have found a way to speed up it by ParallelTable, this is the code:

len = ((Sqrt[1.93*10^16] - 1*10^8)/40 // IntegerPart)*10;
isMatchQ =
If[IntegerDigits[#^2][[1 ;; -1 ;; 2]] ==  {1, 2, 3, 4, 5, 6, 7, 8,
9}, Throw[#]] &;
ParallelTable[n = 10^8 + {3, 7} + len*i;
Catch[Do[isMatchQ /@ n; n += 10, {len / 10}]], {i, 0,
3}] // AbsoluteTiming
(*{42.523432, {Null, Null, Null, 138901917}}*)


I divided the whole range of n which makes 10^16 <= n^2 <= 1.93*10^16 into 4 parts, the length of each part is len. And I calculate these 4 parts in 4 kernels at the same time by ParallelTable, then get the result in a much shorter time, only 42.5s. (Sorry for my poor English....)

Here's a slight improvement of @wuyingddg's Catch&Throw method. The key point is to make Throw working inside ParallelDo using the method mentioned in this post:

SetSharedFunction[pThrow];
pThrow[expr_] := Throw[expr]
isMatchQ =
If[IntegerDigits[#^2][[1 ;; -1 ;; 2]] == {1, 2, 3, 4, 5, 6, 7, 8, 9, 0}, pThrow[#]] &;
Catch[ParallelDo[
isMatchQ /@ (i + {30, 70}), {i, 10^9, Sqrt[193 10^16], 100}]] // AbsoluteTiming
`
• It is much more brief than the ParallelTable one and spend nearly same time, pretty good! – wuyingddg Sep 3 '14 at 14:21