\documentstyle[A4,11pt]{article}
\begin{document}
\title{Computations on the $3n+1$ conjecture}
\author{Gaston H. Gonnet\\Informatik E.T.H. Z\"{u}rich\\Switzerland}
\maketitle
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} = \left\{ \begin{array}{ll}
	3n_i+1 & \mbox{if $n_i \equiv 1 \pmod{2}$}\\ 
	n_i / 2 & \mbox{if $n_i \equiv 0 \pmod{2}$}
	\end{array}
\right.
\]
For example, if $n_0=7$ the sequence starts as
7, 22, 11, 34, 17, 52, 26, 13, 40, 20,
10, 5, 16, 8, 4, 2, 1, 4, 2, 1, 4, ...
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: \verb=DistTo1= computes the number of iterations
that will take its argument to reach the value 1 for the
first time and \verb=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 ) + 2k \]
\[ \sigma ( n 4^k + 1 ) = \sigma ( n 3^k + 1 ) + 3k \]
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 \pmod{2^k}$,
for example,
\[ \sigma ( 8n+3 ) = \sigma( 9n+4 ) + 5 \]
\[ \sigma ( 8n+5 ) = \sigma( 3n+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_{2k}$ and $b_{2k}$ 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_{2k}(n) = \lfloor q/2^k \rfloor 3^{b_k(q \bmod 2^k)} +
a_k(q \bmod 2^k) \]
\[ b_{2k}(n) = b_k(n_0) + b_k(q \bmod 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:
\small\begin{verbatim}
> DistTo1 := proc(n) local its, t, qt, u, rt, bits;
>     if not type([args],[integer]) or n < 1 then
>         ERROR(`invalid arguments`) fi;
>     its := 0;
>     t := n;
>     do
\end{verbatim}
\normalsize
If the argument is small enough (or becomes small enough)
compute the result directly.
\small\begin{verbatim}
>         if t < 10 then RETURN( op(t,[0,1,7,2,5,8,16,3,19,6]) + its ) fi;
\end{verbatim}
\normalsize
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).
\small\begin{verbatim}
>         bits := 3;
>         while 2*length(pow2(bits)) < length(t) do bits := 2*bits od;
\end{verbatim}
\normalsize
The function \verb=iquo= returns the quotient and stores the
remainder in its third argument.
\small\begin{verbatim}
>         qt := iquo(t,pow2(bits),'rt');
>         u := CollatzMod(rt,bits);
\end{verbatim}
\normalsize
The values $a_{bits}$ and $b_{bits}$ are in the list \verb=u=,
now add the number of iterations and compute the next value
of \verb=t=.
\small\begin{verbatim}
>         its := its+u[2]+bits;
>         t := qt*3^u[2]+u[1];
>     od;
> end:
\end{verbatim}
\normalsize

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.
\small\begin{verbatim}
> CollatzMod := proc(r,b)
> local b2, q1, q2, u1, u2, res;
>     b2 := iquo(b,2);
\end{verbatim}
\normalsize
The function \verb=irem= returns the remainder and stores the
quotient in its third argument.
\small\begin{verbatim}
>     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]];
\end{verbatim}
\normalsize
If \verb=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.
\small\begin{verbatim}
>     if b <= 12 then CollatzMod(args) := res fi;
>     res
> end:
\end{verbatim}
\normalsize
Now we will define all the possible cases for 3 or less bits,
this will terminate the recursion of the above function.
(Notice that \verb=b= is always of the form $3 \times 2^k$).
\small\begin{verbatim}
> 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]:
\end{verbatim}
\normalsize

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 \times 2^b}$ it is
worthwhile computing them efficiently and remembering
the past values.
\small\begin{verbatim}
> pow2 := proc(b) option remember;
> if b <= 24 then 2^b else pow2(b/2)^2 fi end:
\end{verbatim}
\normalsize
Now compute some values (and try to verify some of the
above identities)
\small\begin{verbatim}
> DistTo1( 2^200 );
\end{verbatim}
\begin{verbatim}
                                    200

\end{verbatim}
\begin{verbatim}
> DistTo1( 2^100-1 ) = DistTo1( 3^100-1 ) + 200;
\end{verbatim}
\begin{verbatim}
                                1465 = 1465

\end{verbatim}
\begin{verbatim}
> DistTo1( 25*4^100+1 ) = DistTo1( 25*3^100+1 ) + 300;
\end{verbatim}
\begin{verbatim}
                                1562 = 1562

\end{verbatim}
\normalsize

The function \verb=Iterate= will start with a value \verb=x=
perform \verb=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.
\small\begin{verbatim}
> Iterate := proc( x, n )
> local bits, qt, rn, rt, rx, u;
>     if not type([args],[posint,posint]) then
>         ERROR(`invalid arguments`) fi;
>     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');
>              u := CollatzMod(rt,bits);
>              rn := rn-u[2]-bits;
>              rx := qt*3^u[2]+u[1];
>              fi
>         od;
>     rx
> end:
\end{verbatim}
\normalsize
No we can try some examples of this function
\small\begin{verbatim}
> Iterate(12345678901234567890,100);
\end{verbatim}
\begin{verbatim}
                             2790345133793393

\end{verbatim}
\begin{verbatim}
> Iterate( 2^100-1, 2^100 );
\end{verbatim}
\begin{verbatim}
                                     1

\end{verbatim}
\normalsize
The next expression should return 1 for any value of \verb=x=
\small\begin{verbatim}
> x := 7052477465274658974657456924727562;
\end{verbatim}
\begin{verbatim}
                  x := 7052477465274658974657456924727562

\end{verbatim}
\begin{verbatim}
> Iterate( x, DistTo1(x) );
\end{verbatim}
\begin{verbatim}
                                     1

\end{verbatim}
\normalsize
Finally, for some timings, this was run on a DEC 5000 workstation
(notice that $3^{50000}-1$ is an integer 23857 digits long).
\small\begin{verbatim}
> st := time():  DistTo1( 3^50000-1 );  time()-st;
\end{verbatim}
\begin{verbatim}
                                  567858

                                  518.350

\end{verbatim}
\normalsize
\end{document}
