# 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. –  Pickett 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

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 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....) - 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}}*)  adding the final zeroes: (* 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 :) – 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 Most PE problems are designed such that any computing system will buckle under the teraflops usually required to brute force a solution. Not so many flops required for this problem, but some experiments can dramatically reduce the number of tests required. The minimum possible square is 1020304050607080900, and the maximum square is 1929394959697989990, corresponding to a range of integers$n$from 1010101010 to 1389026623. Apparently, there are about 380 million integers to test. However, the only possible square ending in 0 must end in 00. So the last three digits of the square are ...900. Therefore, the integers$n$to be squared must all end in either 30 or 70, and the number of such integers to test is now only about 7.6 million. The required integer is unique, according to the problem statement, so a simple, direct use of ParallelSum with Boole is possible. Avoid MatchQ as @Pickett suggests. AbsoluteTiming[ ParallelSum[ n*Boole[IntegerDigits[n^2][[;; ;; 2]] == {1,2,3,4,5,6,7,8,9,0}], {n, 1010101030, 1390000000, 100}] + ParallelSum[ n*Boole[IntegerDigits[n^2][[;; ;; 2]] == {1,2,3,4,5,6,7,8,9,0}], {n, 1010101070, 1390000000, 100}]] (* {5.915985, 1389019170} *)  Further investigation shows that there are only three possible 5-digit endings for the required square: {{8,0,9,0,0},{8,4,9,0,0},{8,8,9,0,0}}. The corresponding 5-digit integers whose squares have these endings are of the form$i+2500k$, where$i=\{10530,11970,10430,12070,10830,11670\}$, and$k=0,1,2,...35$. The beginning 5 digits of the integer range from$j=10101$to$j=13890$. Thus, integers$n$of the form$n=100000*j+i+2500k$have squares with the required 8_ 9_ 0 ending. Now there are only about 821 thousand integers$n\$ to test. Approximately 460 times fewer tests than a direct brute force approach.

AbsoluteTiming[
ParallelSum[
(100000 j + i + 2500 k) *
Boole[IntegerDigits[(100000j+i+2500k)^2][[;; ;; 2]] == {1,2,3,4,5,6,7,8,9,0}],
{j, 10101, 13890},
{k, 0, 35},
{i, {10530, 11970, 10430, 12070, 10830, 11670}}]]
(* {1.002618, 1389019170} *)


But nothing beats counting down from the maximum as @MarkMcClure comments.

-

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