Computations on the 3 n + 1 conjecture
Worksheet created by Michael Monagan, Text taken from Collatz.tex, Gaston H. Gonnet, Informatik E.T.H. Zurich, Switzerland
Abstract: Worksheet performs number of iterations in Collatz's problem, i.e., the 3x+1 conjecture.
Application Areas/Subjects: Number Theory
Keywords: Collatz, number Theory, Ulam, Syracuse, Kakutani, Hasse
> restart:
The 3n+1 sequence has probably consumed more CPU time than any other number theoretic conjecture. The reason being that its statement is so simple, that most amateurs will feel compelled to write programs to test it. This sequence, attributed to Lothar Collatz, has bee given various names, including Ulam's conjecture, Syracuse problem, Kakutani's problem and Hasse's algorithm. The conjecture is based on the iteration defined by:
n[i+1] = if n[i] is even then n[i]/2 else 3*n[i]+1 .
For example, if n[1]=7 the sequence goes as follows
> n[1] := 7;
> for i from 1 to 20 do
> if irem(n[i],2)=0 then n[i+1] := iquo(n[i],2) else n[i+1] := 3*n[i]+1 fi;
> od;
> seq( n[i], i=1..20 );
In this case it is obvious that the sequence entered a cycle and will repeat the values 4,2,1 forever. The 3n+1 conjecture states that for any starting positive integer, the sequence eventually reaches 1 and repeats as above.
Jeffrey Lagarias published an excellent article summarizing most of the work done on this problem. His article appeared in the American Mathematical Monthly, Volume 92, Number 1, January 1985, pages 3-23. Any reader who is interested in doing work on this problem is urged to read Lagarias' article.
The purpose of this note is to show how to compute iterations of this series efficiently using Maple. The programs that we will show, can compute thousands of iterations over terms which are thousands of digits long in reasonably short times.
We provide two main functions for working with these sequences: ``DistTo1'' computes the number of iterations that will take its argument to reach the value 1 for the first time and ``Iterate'' iterates its first argument the number of times specified in its second argument. Both functions use a recursive modular computation of the terms which requires some additional explanation.
Let sigma(n) denote the number of iterations required to reduce n to 1, or what we call the distance to 1.
The following relations are known and/or easy to derive:
> sigma( n*2^k ) = sigma (n) + k;
> sigma ( n*2^k - 1 ) = sigma( n*3^k - 1 ) + 2*k;
> sigma ( n*4^k + 1 ) = sigma ( n*3^k + 1 ) + 3*k;
It is easy to see that for any value of n, the first k iterations which divide by 2 depend exclusively on the value of n mod 2^k, for example.
> sigma( 8*n+3 ) = sigma( 9*n+4 ) + 5;
> sigma( 8*n+5 ) = sigma( 3*n+2 ) + 4;
It should also be observed that the factor that multiplies n in the resulting term is always a power of 3. Furthermore, this power of 3 corresponds exactly to the number of 3n+1 iterations necessary while doing the k divide by two iterations.
The above suggests that we could run the iteration in symbolic terms, say for an arbitrary value of k by computing the coefficients for the resulting iterations. This strategy has two immediate advantages. First, for very large numbers, we will need to do the computation with long numbers only once for every k division by two steps; secondly, we could easily tabulate all the possible values for small values of k, and hence save in the computation of the coefficients themselves. In mathematical terms, what we will compute are the functions a_k and b_k defined by the following equation
> sigma ( 2^k*n[1] + n[0] ) = sigma ( 3^(b[k]( n[0] ))*n[1] +
> a[k]( n[0] ) ) + k + b[k]( n[0] );
By analyzing two steps of such an iteration we can derive the function a[2*k] and b[2*k] in terms of a[k] and b[k]. Let
> n = n[1]*2^k + n[0];
> q = n[1]*3^(b[k](n[0])) + a[k](n[0]);
then
> a[2*k](n) = floor( q/2^k ) * 3^(b[k](' q mod 2^k ')) + a[k]('q mod 2^k');
> b[2*k](n) = b[k](n[0]) + b[k]('q mod 2^k');
As it turns out, this gives us a very efficient method for computing several iterations of the recurrence at once. The main program to compute the sigma(n) function will compute k as an appropriate power of two, as large as possible, but guaranteeing that the series will not reach 1 in less than k steps. This is achieved simply by insuring that 2^k < n. The function definition, argument checking and initialization look like:
> DistTo1 := proc(n::posint) local its, t, qt, u, rt, bits;
> its := 0;
> t := n;
> do
If the argument is small enough (or becomes small enough) compute the result directly.
> if t < 10 then RETURN( op(t,[0,1,7,2,5,8,16,3,19,6]) + its ) fi;
Now compute 2^bits < t by comparing the lengths of their decimal representations (the iteration is stopped before exceeding t so that we save one (the largest) power of two computation).
> bits := 3;
> while 2*length(pow2(bits)) < length(t) do bits := 2*bits od;
The function iquo returns the quotient and stores the remainder in its third argument.
> qt := iquo(t,pow2(bits),'rt');
> u := CollatzMod(rt,bits);
The values a[bits] and b[bits] are in the list u, now add the number of iterations and compute the next value of t.
> its := its+u[2]+bits;
> t := qt*3^u[2]+u[1];
> end:
The function which computes a[bits] and b[bits] is written recursively. The bottom of the recursion is entered as fixed constants at the end of the definition. The body of this function just follows the above mathematical definition.
> CollatzMod := proc(r,b)
> local b2, q1, q2, u1, u2, res;
> b2 := iquo(b,2);
The function irem returns the remainder and stores the quotient in its third argument.
> u1 := CollatzMod(irem(r,pow2(b2),'q1'),b2);
> u2 := CollatzMod(irem(q1*3^u1[2]+u1[1],pow2(b2),'q2'),b2);
> res := [q2*3^u2[2]+u2[1], u1[2]+u2[2]];
If b is less or equal to 12, then we will remember any computed value for the future. This will effectively create a dynamic table (will only grow as needed) with at most 4175=2^12+2^6+2^3+2^2+2^1 entries.
> if b <= 12 then CollatzMod(args) := res fi;
> res;
Now we will define all the possible cases for 3 or less bits, this will terminate the recursion of the above function. (Notice that b is always of the form 3 2^k).
> CollatzMod(0,1) := [0,0]: CollatzMod(1,1) := [2,1]:
> CollatzMod(0,2) := [0,0]: CollatzMod(1,2) := [1,1]:
> CollatzMod(2,2) := [2,1]: CollatzMod(3,2) := [8,2]:
> CollatzMod(0,3) := [0,0]: CollatzMod(1,3) := [2,2]:
> CollatzMod(2,3) := [1,1]: CollatzMod(3,3) := [4,2]:
> CollatzMod(4,3) := [2,1]: CollatzMod(5,3) := [2,1]:
> CollatzMod(6,3) := [8,2]: CollatzMod(7,3) := [26,3]:
The above functions will use some powers of two again and again. Since these powers of two are very few (just the ones of the form 2^(3 2^b) it is worthwhile computing them efficiently and remembering the past values.
> pow2 := proc(b) option remember;
> if b <= 24 then 2^b else pow2(b/2)^2 fi end:
Now compute some values (and try to verify some of the above identities)
> DistTo1( 2^200 );
> DistTo1( 2^100-1 ) = DistTo1( 3^100-1 ) + 200;
> DistTo1( 25*4^100+1 ) = DistTo1( 25*3^100+1 ) + 300;
The function Iterate will start with a value x perform n Collatz iterations and return the resulting term. It is written using the same ideas of modular computation as above, plus a provision for terminating quickly once that the iteration reaches its cycle at 4,2,1.
> Iterate := proc( x::posint, n::posint )
> local bits, qt, rn, rt, rx, u;
> rn := n; rx := x;
> while rn > 0 do
> if rx=1 then RETURN( op( modp(rn,3)+1, [1,4,2] ) )
> elif rx=2 then RETURN( op( modp(rn,3)+1, [2,1,4] ) )
> elif rx=4 then RETURN( op( modp(rn,3)+1, [4,2,1] ) )
> elif rn < 6 then
> if irem(rx,2,'qt')=1 then rx := 3*rx+1 else rx := qt fi;
> rn := rn-1
> else bits := 3;
> while 2*length(pow2(bits)) < length(t) and 4*bits < rn
> do bits := 2*bits od;
> qt := iquo(rx,pow2(bits),'rt');
> rn := rn-u[2]-bits;
> rx := qt*3^u[2]+u[1];
> fi
> rx;
>
Now we can try some examples of this function
> Iterate(12345678901234567890,100);
> Iterate( 2^100-1, 2^100 );
The next expression should return 1 for any value of x
> x := 7052477465274658974657456924727562;
> Iterate( x, DistTo1(x) );
Finally, for some timings, this was run on a DEC 5000 workstation (notice that 3^50000-1 is an integer 23857 digits long).
> st := time(): DistTo1( 3^5000-1 ); time()-st;
> st := time(): DistTo1( 3^50000-1 ); time()-st;
