# 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

• Take a look at memoization Oct 3, 2013 at 15:01

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


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

• What's the usage of s ? Sep 20, 2014 at 7:06
• @xzczd I forgot to include the plotting command. Edited. Sep 20, 2014 at 20:27
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

• Wow. Recursive functions always make me feel like aliens have landed. Jan 22, 2015 at 14:48