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

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

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"];