Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
share|improve this question
    
Take a look at memoization –  belisarius Oct 3 '13 at 15:01

2 Answers 2

up vote 11 down vote accepted

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"];
share|improve this answer

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}

share|improve this answer
    
What's the usage of s ? –  xzczd yesterday
    
@xzczd I forgot to include the plotting command. Edited. –  belisarius 15 hours ago

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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