4
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ Take a look at memoization $\endgroup$ Oct 3, 2013 at 15:01

3 Answers 3

11
$\begingroup$

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"];
$\endgroup$
4
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ What's the usage of s ? $\endgroup$
    – xzczd
    Sep 20, 2014 at 7:06
  • 1
    $\begingroup$ @xzczd I forgot to include the plotting command. Edited. $\endgroup$ Sep 20, 2014 at 20:27
2
$\begingroup$
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
$\endgroup$
1
  • $\begingroup$ Wow. Recursive functions always make me feel like aliens have landed. $\endgroup$
    – yippy_yay
    Jan 22, 2015 at 14:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.