Calculating Collatz Cousins

Question

Define the function f(n) for a positive integer n as follows:

• n / 2, if n is even
• 3 * n + 1, if n is odd

If you repeatedly apply this function to any n greater than 0, the result always seems to converge to 1 (though nobody's been able to prove that yet). This property is known as the Collatz Conjecture.

Define an integer's stopping time as the number of times you have to pass it through the Collatz function f before it reaches 1. Here are the stopping times of the first 15 integers:

1  0
2  1
3  7
4  2
5  5
6  8
7  16
8  3
9  19
10 6
11 14
12 9
13 9
14 17
15 17

Let's call any set of numbers with the same stopping time Collatz cousins. For example, 5 and 32 are Collatz cousins, with a stopping time of 5.

Your task: write a program or function that takes a nonnegative integer and generates the set of Collatz cousins whose stopping time is equal to that integer.

Input

A nonnegative integer S, given via STDIN, ARGV, or function argument.

Output

A list of all numbers whose stopping time is S, sorted in ascending order. The list may be output by your program, or returned or output by your function. Output format is flexible: space-separated, newline-separated, or any standard list format of your language is fine, as long as the numbers are easily distinguishable from one another.

Requirements

Your submission must give correct results for any S ≤ 30. It should finish in seconds or minutes, not hours or days.

Examples

0  -> 1
1  -> 2
5  -> 5, 32
9  -> 12, 13, 80, 84, 85, 512
15 -> 22, 23, 136, 138, 140, 141, 150, 151, 768, 832, 848, 852, 853, 904, 906, 908, 909, 5120, 5376, 5440, 5456, 5460, 5461, 32768

Here is a Gist of the output for S = 30.

This is code-golf: shortest program in bytes wins. Good luck!

Answer 1

Mathematica, 98 92 89 bytes

This solution solves S = 30 immediately:

(p={0};l={1};Do[l=Complement[##&@@{2#,Mod[a=#-1,2]#~Mod~3~Mod~2a/3}&/@l,p=p⋃l],{#}];l)&

This is an unnamed function taking S as its only parameter and returning a list of the Collatz cousins.

The algorithm is a simple breadth-first search. The Collatz cousins for a given S are all the integers that can be reached from the Collatz cousins for S-1 via 2*n or odd numbers that can be reached via (n-1)/3. We also need to ensure that we only produce those integers which were reached for the first time after S steps, so we keep track of all previous cousins in p and remove those from the result. Since we're doing that anyway, we can save a few bytes by computing the steps from all previous cousins (not just those from S-1) to save a few bytes (that makes it slightly slower, but not noticeably for the required S).

Here is a slightly more readable version:

(
  p = {0};
  l = {1};
  Do[
    l = Complement[
      ## & @@ {2 #, Mod[a = # - 1, 2] #~Mod~3~Mod~2 a/3} & /@ l,
      p = p ⋃ l
    ]~Cases~_Integer,
    {#}
  ];
  l
) &

Answer 2

Pyth, 26 24 21 bytes

Su+yMG-/R3fq4%T6G1Q]1

This code runs instantly for S=30. Try it out yourself: Demonstration

Thanks to @isaacg for saving 5 bytes.

Explanation

My code starts with 1 and undos the Collatz function. It maps all numbers d of the S-1 step to 2*d and (d-1)/3. The last one in not always valid though.

                        implicit: Q = input number
                   ]1   start with G = [1]
 u                Q     apply the following function Q-times to G:
                          update G by
   yMG                      each number in G doubled
  +                       +
          fq4%T6G           filter G for numbers T, which satisfy 4==T%6
       /R3                  and divide them by 3
      -          1          and remove 1, if it is in the list
                            (to avoid jumping from 4 to 1)
S                       sort the result and print

Answer 3

Python 2, 86 83 75 73 71 bytes

f=lambda n,k=1:sorted([k][n:]or(k>4==k%6and f(n-1,k/3)or[])+f(n-1,k*2))

Call like f(30). n = 30 is pretty much instant.

(Thanks to @DLosc for the idea of recursing by k being a number rather than a list of cousins, and a few bytes. Thank to @isaacg for dropping ~-.)

This variant is much shorter, but unfortunately takes too long due to exponential branching:

f=lambda n,k=1:sorted([k][n:]or(k>4==k%6)*f(n-1,k/3)+f(n-1,k*2))

Answer 4

CJam, 29 26 bytes

Xari{{2*_Cmd8=*2*)}%1-}*$p

Credit goes to @isaacg for his idea to remove 1's after each iteration, which saved me two bytes directly and another one indirectly.

Try it online in the CJam interpreter (should finish in less than a second).

How it works

Xa       e# Push A := [1].
ri{      e# Read an integer from STDIN and do the following that many times:
  {      e# For each N in A:
    2*   e#     Push I := (N * 2) twice.
    _Cmd e#     Push (I / 12) and (I % 12).
     8=  e#     Push K := (I % 12 == 8).

         e#     (K == 1) if and only if the division ((N - 1) / 3) is exact and
         e#     yields an odd integer. In this case we can compute the quotient 
         e#     as (I / 12) * 2 + 1.

    *2*) e#     Push J := (I / 12) * K * 2 + 1.

         e#     This yields ((N - 1) / 3) when appropriate and 1 otherwise.
  }%     e# Replace N with I and J.
  1-     e# Remove all 1's from A.

         e# This serves three purposes:

         e# 1. Ones have been added as dummy values for inappropriate quotients.

         e# 2. Not allowing 1's in A avoids integers that have already stopped
         e#    from beginning a new cycle. Since looping around has been prevented,
         e#    A now contains all integers of a fixed stopping time.

         e# 3. If A does not contain duplicates, since the maps N -> I and N -> J
         e#      are inyective (exluding image 1) and yield integers of different
         e#      parities, the updated A won't contain duplicates either.

}*       e#
$p       e# print(sort(C))

Answer 5

CJam, 35 bytes

1]ri{_"(Z/Y*"3/m*:s:~\L+:L-_&0-}*$p

Explanation coming soon. This is a much faster version than the "pretty straight forward" approach (see it in edit history).

Try it online here for N = 30 which runs in seconds on the online version and instantly in the Java Compiler

Answer 6

Java 8, 123

x->java.util.stream.LongStream.range(1,(1<<x)+1).filter(i->{int n=0;for(;i>1;n++)i=i%2<1?i/2:3*i+1;return n==x;}).toArray()

When x is 30, the program takes 15 minutes and 29 seconds.

Expanded

class Collatz {
    static IntFunction<long[]> f =
            x -> java.util.stream.LongStream.range(1, (1 << x) + 1).filter(i -> {
                int n = 0;
                for (; i > 1; n++)
                    i = i % 2 < 1 ? i / 2 : 3 * i + 1;
                return n == x;
            }).toArray();

    public static void main(String[] args) {
        System.out.println(Arrays.toString(f.apply(15)));
    }
}

Answer 7

Python 2, 118 bytes

Well, I figured that I wouldn't reach the best Python score after seeing @Sp3000's solution. But it looked like a fun little problem, so I wanted to try an independent solution anyway:

s={1}
for k in range(input()):
 p,s=s,set()
 for t in p:s.add(2*t);t>4and(t-1)%6==3and s.add((t-1)/3)
print sorted(s)

Same thing before stripping whitespace:

s={1}
for k in range(input()):
    p,s=s,set()
    for t in p:
        s.add(2 * t)
        t > 4 and (t - 1) % 6 == 3 and s.add((t - 1) / 3)
print sorted(s)

This is a very direct implementation of a breadth first search. In each step, we have the set with stopping time k, and derive the set with stopping time k + 1 by adding the possible predecessors of each value t in the set from step k:

• 2 * t is always a possible predecessor.
• If t can be written as 3 * u + 1, where u is an odd number that is not 1, then u is a predecessor as well.

Takes about 0.02 seconds to run for N = 30 on my MacBook Pro.

Answer 8

PHP 5.4+, 178 bytes

The function

function c($s,$v=1,$p=[],&$r=[]){$p[]=$v;if(!$s--){return$r[$v][]=$p;}c($s,$v*2,$p,$r);is_int($b=($v-1)/3)&!in_array($b,$p)&$b%2?c($s,$b,$p,$r):0;ksort($r);return array_keys($r);}

Test & Output

echo "0 - ".implode(',',c(0)).PHP_EOL;
// 0 - 1
echo "1 - ".implode(',',c(1)).PHP_EOL;
// 1 - 2
echo "5 - ".implode(',',c(5)).PHP_EOL;
// 5 - 5,32
echo "9 - ".implode(',',c(9)).PHP_EOL;
// 9 - 12,13,80,84,85,512
echo "15 - ".implode(',',c(15)).PHP_EOL;
// 15 - 22,23,136,138,140,141,150,151,768,832,848,852,853,904,906,908,909,5120,5376,5440,5456,5460,5461,32768

S(30) runs in 0.24 seconds*, returns 732 elements. A couple are

86,87,89,520,522,524,525,528, [ ... ] ,178956928,178956960,178956968,178956970,1073741824

*To save on bytes, I had to add ksort and array_keys at every step. The only other choice I had was to make a small wrapper function that calls c() and then calls array_keys and ksort on the result once. But due to the time still being decently snappy, I decided to take the performance hit for low byte count. Without the proper sorting & processing, the time is 0.07 seconds on average for S(30).

If anyone has any clever ways of getting the proper processing only once without too many additional bytes, please let me know! (I store my numbers as array keys, hence the use of array_keys and ksort)

Answer 9

tinylisp, 993 823 bytes

Shout out to @DLosc for carving off 170 bytes! Woo Hoo!

(d !(q(1 1
(d B(q(S S
(d C(q((S T)(i T(C(c(a(h S)(h S))S)(s T 1))S
(d F(q((S T)(i T(F(c((q((U)(a(a U U)U)))(h(t S)))S)(s T 1))S
(d G(F ! 38
(d H(q((S T)(i(e T !)S(i(l(h T)(a S 1))(H(s S(h T))(t T))(H S(t T
(d K(q((S T)(i(e T !)S(i(l(h T)(a S 1))(a(h(t(t T)))(K(s S(h T))(t T)))(K S(t T
(d M(q((S)(N(i(a(H(s S 1)G)(a(e S 4)(H S(C ! 30))))()(B(K(s S 1)G)))(B(a S S
(d N(q((S T)(i S(N(t S)(c(h S)T))T
(d O(q((T S)(i S(O(N(M(h S))T)(t S))T
(d P(q((T S)(i S(Q(c(h S)T)(t S))T
(d Q(q((T S)(i S(P T(t S))T
(d R(q((S T U)(i S(i T(i(e(h S)(h T))(R(t S)(t T)(c(h S)U))(i(l(h S)(h T))(R(t S)T(c(h S)U))(R S(t T)(c(h T)U))))(R(t S)T(c(h S)U)))(i T(R S(t T)(c(h T)U))U
(d V(q((S T)(i S(V(t S)(c(h S)T))T
(d W(q((S)(i(t S)(V(R(W(Q()S))(W(P()S))())())S
(d Y(q((T S)(i S(Y(i T(O()T)())(s S 1))(B T S
(d A(q((S)(i S(W(h(Y(B 2)(s S 1))))(B 1

Try it online!

Ungolfed:

(d list
(q(args
  args
)))


(d make-pow2
(q((S T)
  (i
    T
    (make-pow2 (c (a (h S) (h S)) S) (s T 1))
    S
  )
)))

(d pow2 (make-pow2 (q(1 1)) 30))

(d times3
(q((U )
   (a (a U U) U)
)))

(d make-pow3
(q((S T)
  (i
    T
    (make-pow3 (c (times3 (h (t S))) S) (s T 1))
    S
  )
)))

(d pow3 (make-pow3 (q (1 1)) 38))

(d trim
(q((S T)
  (i
    (e T (q(1 1)))
    S  
    (i 
      (l (h T) (a S 1))
      (trim (s S (h T)) (t T))
      (trim S (t T))
    )
  )
)))

(d div3?
(q((S)
  (e 0 (trim S pow3))
)))

(d div2?
(q((S)
  (e 0 (trim S pow2))
)))


(d div3*
(q((S T)
  (i
    (e T (q(1 1)))
    S
    (i 
      (l (h T) (a S 1))
      (a (h (t (t T))) (div3* (s S (h T)) (t T)))       
      (div3* S (t T)) 
    )
  )  
)))

(d div3
(q((S)
  (div3* S pow3)
)))




(d G
(q((S)
  (i 
    (div3? (s S 1))
    (i 
      (e (div3 (s S 1)) 1)
      (list (a S S))
      (i
        (div2? (div3 (s S 1)))
        (list (a S S))
        (list (div3 (s S 1)) (a S S))
      )      
    )
    (list (a S S))
  )
)))

(d concat
(q((S T)
  (i
    S
    (concat (t S) (c (h S) T))
    T
  )
)))

(d step*
(q((N acc)
  (i
    N
    (step* (t N) (concat (G (h N)) acc )) 
    acc
  ) 
)))

(d even
(q((S T)
  (i
    (e S ())
    T
    (odd (t S) (c (h S) T))        
  )    
)))

(d odd
(q((S T)
  (i
    (e S ())
    T
    (even (t S) T)        
  )    
)))

(d merge
(q((S T U)
  (i
    (e S ())
    (i
      (e T ())
      U
      (merge S (t T) (c (h T) U)) 
    )
    (i
      (e T ())
      (merge (t S) T (c (h S) U)) 
      (i
        (e (h S) (h T))
        (merge (t S) (t T) (c (h S) U))      
        (i
          (l (h S) (h T))
          (merge (t S) T (c (h S) U))
          (merge S (t T) (c (h T) U))
        )
      )
    )
  )  
)))

(d rev
(q((S T)
  (i 
    (e S ())
    T
    (rev (t S) (c (h S) T))
  )
)))

(d sort
(q((S)
  (i
    (e () S)
    S
    (i 
      (e () (t S))
      S
      (rev (merge (sort(odd S ())) (sort (even S ())) ()) ())
    )
  )
)))

(d step
(q((N)
  (i
    N
    (step* N ())  
    ()  
  )
)))

(d stopping
(q((N S)
  (i
    (e S 0)
    (list N S)
    (stopping (step N) (s S 1))
  )  
)))

(d A
(q((S)
  (i
    (e S 0)
    (q(1))
    (sort (h (stopping (q(2)) (s S 1))) )
  )  
)))

Try it online!

Explanation

My initial approach was a naive one where I started from 2^n and counted down checking each number to see if it's stopping number was n. That was way too slow. After consulting other answers for inspiration, I switched to using the reverse Collatz function, generating the tree, and stopping after n iterations.

I discovered by experimenting, that the tree does not grow exponentially - although n*2 is always a hit, the (n-1)/3 path only hits occasionally. So by the time we get to n=30 we have about twice the number of integers that we want, but no more than that.

As a result, I put the sort and de-dup at the very end, and called it just once - this saves some execution time.

The sort is a merge sort - this makes it convenient to de-dup while merging.

The arithmetic operations required: divides-by-two? divides-by-three? and divide-by-three all had to be optimized - the naive implementations of these were too slow. To optimize, I created a list of powers of 2 and 3, and used these to make the operations faster.

Testing

My code has been tested and verified for n = 0,1,5,9,15 and 30. All are correct. n=30 runs in about 15 seconds on TIO.

Answer 10

C Language

#include <stdio.h>
#include <limits.h>    
const int s = 30;

bool f(long i)
{
    int r = 0;
    for(;;)
        if (i < 0 || r > s) return false;
        else if (i == 1) break;
        else{r ++;i = i % 2 ? 3*i + 1 : i/2;}
    return (r==s);
}

void main(){
    for(long i = 1; i < LONG_MAX; i++) if (f(i)) printf("%ld ", i);
}