How to speed up my Project Euler code

Question

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.

Answer 1

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

Answer 2

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

Answer 3

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

Answer 4

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