This is a continuation of my previous post on Collatz Lengths.
Let me start by repeating the function that I would like to change a bit. Originally I had \[ f(n) = \begin{cases} n/2 & \text{ if $n$ is even }\\ 1 & \text{ if $n=1$ }\\ 3n+1 & \text{ if $n\neq 1$ is odd } \end{cases} \] I used this function to define a recursive sequence \(a_{n+1} = f(a_n)\). And the question was to show that this sequence eventually produces \(1\) for any initial natural number value \(a_0\).
For a prime number \(p\) and a natural number \(n\), the \(p\)-part of \(n\) is \(\nu_p(n)\) the largest natural number \(m\) such that \(p^m\) divides \(n\).
Now, notice that the recursive sequence gradually destroys the 2-part of any initial value first, then applies the function \(3x+1\), and it repeats the procedure until it hits \(1\). This is where we stop.
Let me concentrate on the crucial part then. I will redefine \(f(n)\) as follows: \[ f(n) = \begin{cases} f(n/2) & \text{ if $n$ is even }\\ 1 & \text{ if $n=1$ }\\ 3n+1 & \text{ if $n\neq 1$ and it is odd } \end{cases} \] Again define a recursive sequence \(a_{n+1} = f(a_n)\), and let \(iter(a)\) be the length of the recursive sequence \((a_n)\) for \(a_0 = a\).
Today, I would like to look the distribution of \(iter(n)/log(n)\).
First, our basic function
(defun f(n)
(cond ((evenp n) (f (/ n 2)))
((= n 1) 1)
(t (1+ (* 3 n)))))
F
Our iteration:
(defvar upper-limit (expt 2 22))
(defun iterate (n &optional (k 0))
(if (= n 1)
k
(iterate (f n) (1+ k))))
ITERATE
OK. Let me run the code.
(with-open-file (results "collatz3.csv"
:direction :output
:if-exists :supersede
:if-does-not-exist :create)
(loop for i from 2 to upper-limit
do (format results "~A ~A~%" i (/ (iterate i) (log i)))))
NIL
Here is the plot of the results (\(n\) versus \(len(n)/log(n)\)):
It is difficult to understand what is going on here. Let me consider the first 20000 values only:
And now for the last 20000 values:
X <- read.csv("collatz3.csv",sep=" ",header=FALSE)
png(filename="collatz3-histogram.png")
hist(X$V2,main="Length(a)/log(a)",xlab="Value")
dev.off()
The distribution of the results is as follows: