The hailstone sequence is a sequence of positive integers from an initial value that is popular through the Collatz conjecture, an unsolved problem in mathematics. We are interested in simulating the sequence at various initial values to gain insights to how the stopping time property relates to the initial value.

Although this property has been studied with respect to the Collatz conjecture, we are primarily interested in the exercise of visualizing the structural properties of sequences using tools provided by the Julia language. The stopping time property has no known analytic solution, so numerical simulation is a method that offers insight.

We begin this experiment with a precise definition of the hailstone sequence and its stopping time property. Then, we proceed to define an appropriate implementation using Julia. Finally, we plot our results using the Gadfly library.

Hailstone Sequence

The hailstone sequence begins with an initial value $a_0$ and subsequent values are determined by

$$ a_i = \begin{cases} \frac{a_{i-1}}{2}, &\text{if $a_{i-1} \equiv 0 \pmod{2}$} \\ 3 a_{i-1} + 1, &\text{otherwise.} \\ \end{cases} $$

Stopping Time

The stopping time is defined as the number of iterations until the number $1$ is reached. For example, if we began the sequence at $a_0 = 4$, then the following sequence values would be

$$ \begin{array}{l} a_1 = 2 \\ a_2 = 1 \end{array} $$

and the stopping time would be $2$. Since the stopping time has no known closed form, we gain some intuition regarding its relation to the initial values by numerical simulation.

Simulating with Julia

We are interested in how the stopping times relate to the initial values. We analyze this relationship by simulating the sequences in Julia and measuring the stopping times explicitly for a small subset of the data $x \in [1, 10000]$.

using Gadfly

# Define the sequence rule
next_state(a) = a % 2 == 0 ? div(a, 2) : 3a + 1

# Define the stopping time
stopping_time(a) = a == 1 ? 0 : 1 + stopping_time(next_state(a))

# Perform 10,000 simulations
xs = 1:10_000
ys = [stopping_time(x) for x in xs]

# Plot the simulation results
plot(x=xs, 
     y=ys, 
     Guide.xlabel("Initial Value"), 
     Guide.ylabel("Stopping Time"))

# Plot the density 
plot(x=ys, 
     Geom.density,
     Guide.xlabel("Stopping Time"))

This code simulates $10,000$ sequences within the range $[1, 10000]$ and records the observed stopping times. The first graph observes the relationship between the stopping time and the initial value directly. The overall graph is increasing at a slow rate.

Stopping Time

The second graph is a plot of the density of the stopping times. The frequency of large stopping times appear to decrease as the value grows larger.

Stopping Time Histogram

Conclusion

Julia provides several tools for rapidly prototyping similar experiments and visualizing the results. In this particular scenario where an analytic solution did not exist for a property of the sequence, numerical simulations demonstrate a capacity for providing insight.