Improve the performance of solutions to Project Euler (#14)

Question

Problem

I have following code, when n<=10^5 it's not slow, but n>2*10^5 it's became very slow. I think maybe some temp value greater than 2^31-1, so make compile invalid. Could you give any idea make it can be compiled? As far as possible not to use recursive algorithm.

pe14 = Compile[{},
   Module[{n1, len, maxLen = 0, res = 0},
    Do[n1 = n;
     len = 1;
     While[n1 != 1,
      n1 = If[EvenQ@n1,  n1~Quotient~2, 3 n1 + 1];
      len++
      ];
     If[len > maxLen, maxLen = len; res = n],
     {n, 1, 10^6}];
    {maxLen, res}
    ]
   ];

pe14[] // AbsoluteTiming

Answer 1

A print statement shows that this will overflow on platforms where Mathematica machine integers are 32 bits.

pe14 = Compile[{}, Module[{n1, len, maxLen = 0, res = 0, print = 0},
    Do[n1 = n;
     len = 1;
     While[n1 != 1, n1 = If[EvenQ@n1, n1~Quotient~2, 3 n1 + 1];
      If[n1 > 10^4*n && print < 10, print++; Print[{n, n1}]];
      len++];
     If[len > maxLen, maxLen = len; res = n], {n, 1, 2 10^5}];
    {maxLen, res}]];

pe14[] // AbsoluteTiming

During evaluation of In[356]:= {77671,1047216490}

During evaluation of In[356]:= {77671,1570824736}

During evaluation of In[356]:= {77671,785412368}

During evaluation of In[356]:= {103561,1047216490}

During evaluation of In[356]:= {103561,1570824736}

During evaluation of In[356]:= {113383,1654740898}

During evaluation of In[356]:= {113383,2482111348}

During evaluation of In[356]:= {113383,1241055674}

During evaluation of In[356]:= {113383,1861583512}

During evaluation of In[356]:= {113383,1325287492}

Out[357]= {4.547872, {383, 156159}}

On 64 bit platforms it will run to completion. Takes around 18 seconds on my desktop.

You can cut a factor of 2 by only explicitly handling odd values, adjusting for evens that are multiples thereof. Another factor of 4 comes from compiling to C.

pe14b = Compile[{{top, _Integer}}, 
   Module[{n1, len, maxLen = 0, res = 0},
    Do[n1 = n;
     len = 1;
     While[n1 != 1, n1 = If[EvenQ@n1, n1~Quotient~2, 3 n1 + 1];
      len++];
     If[len > maxLen, maxLen = len + Floor[Log[2, top/N[n]]]; 
      res = n], {n, 1, top, 2}];
    {maxLen, res}], CompilationTarget -> "C"];

In[377]:= pe14b[10^6] // AbsoluteTiming

Out[377]= {2.396513, {525, 837799}}

For platforms that do not support 64 bit machine integers, one might try to emulate this with machine doubles. EvenQ would have to check whether dividing by two creates a fractional part. The variant below seems to work as it ought.

pe14c = Compile[{{top, _Integer}}, 
   Module[{n1, len, maxLen = 0, res = 0.},
    Do[n1 = n;
     len = 1;
     While[n1 != 1, 
      n1 = If[FractionalPart[n1/2.] == 0, n1/2, 3 n1 + 1];
      len++];
     If[len > maxLen, maxLen = len + Floor[Log[2, top/n]]; 
      res = n], {n, 1., top, 2.}];
    {maxLen, res}], CompilationTarget -> "C"];

Answer 2

Here is a (tolerably slow) non compiled version I wrote some years ago:

s = Array[{}, 10^6];
f[x_] := f[x] = If[x == 1, 1, 1 + f[If[EvenQ[x], x/2, (3 x + 1)]]]; 
Block[{$RecursionLimit = 1000, a = 0, j},
     Do[If[a < (s[[i]] = f[i]), a = f[i]; j = i], {i, Reverse@Range@1000000}];
       Print@a; Print@j;
       ]

{Max@s, ListPlot@Tally@s}

{525, }

Answer 3

ClearAll[col]
col[1]=0
col[n_] := col[n] = If[EvenQ[n], col[n/2], col[3 n + 1]] + 1
ctab = Table[col[n], {n, 1, 1000000}]; // AbsoluteTiming
  {13.358787, Null}
Max[ctab]
  524