Sequential Dependence and Residual Structure in Collatz Stopping Times: Empirical Tests of the Stochastic Model

Abstract

The standard stochastic model for the Collatz map, originating with Terras (1976) and developed by Lagarias and Kontorovich (2009), predicts that stopping times grow as approximately 6.95 log n with parity sequences behaving like independent coin flips. We test this model with three quantitative experiments on the first 10^7 positive integers. First, we show that the mean stopping time increment per trailing 1-bit in binary is 6.17 steps, matching the heuristic prediction of 6.23 steps to within 1.0%, providing a precise quantitative validation of the probabilistic framework. Second, we measure the autocorrelation function of stopping times at lags from 1 to 1024 and find that it follows a power law r(h) ~ 0.321 h^{-0.141} (R^2 = 0.884), decaying far more slowly than any exponential model can accommodate (exponential R^2 = 0.315). This reveals long-range sequential dependence in the stopping-time landscape, arising from the exact identity st(2m) = 1 + st(m). Third, we test conditional independence by computing within-class autocorrelation for residues modulo 2^k and find that each conditioning bit reduces the within-class lag-1 autocorrelation by a nearly constant decrement, from 0.396 at k = 1 to 0.106 at k = 6, without reaching zero. The binary suffix captures variance at a rate of exactly 1.00% per bit, yet the autocorrelation structure is not fully removable by finite conditioning. These results quantify the gap between the independent-parity idealization and the deterministic arithmetic that generates it.

1. Introduction

The Collatz function C: N -> N maps even integers to n/2 and odd integers to 3n + 1. The total stopping time sigma(n) counts the number of iterations required for n to reach 1. The Collatz conjecture asserts that sigma(n) is finite for every n >= 1. While the conjecture remains unresolved, the statistical behavior of stopping times is well described by stochastic models that treat the parity of successive iterates as independent fair coin flips (Terras, 1976; Lagarias, 1985; Kontorovich and Lagarias, 2009; Tao, 2019).

In the standard heuristic, each iteration replaces n by approximately (3/4)n on average, predicting a mean total stopping time of log(n)/log(4/3) ~ 6.95 log n. This prediction is remarkably accurate for large n. However, individual stopping times are deterministic functions of n, and adjacent integers can share trajectory fragments, creating dependencies that the independent-parity model ignores.

This paper investigates the fine structure of these dependencies through three complementary experiments. Section 2 establishes notation and reviews the probabilistic framework. Section 3 validates the stochastic model by computing the per-bit stopping time increment from trailing binary structure. Section 4 measures the autocorrelation function of the stopping-time sequence and identifies its decay law. Section 5 tests whether conditioning on the binary suffix removes the sequential dependence. Section 6 reports supplementary results on the growth of maximal stopping times. Section 7 discusses implications for the stochastic model and open questions.

All computations use the first N = 10^7 positive integers with Python's standard library and random.seed(42) for reproducibility.

2. Notation and the Stochastic Framework

Define the total stopping time sigma(n) = min{k >= 0 : C^k(n) = 1}, where C^k denotes k-fold iteration. For n >= 2, sigma(n) >= 1. Let v_2(n) denote the 2-adic valuation of n, the largest power of 2 dividing n, and let tau(n) denote the number of trailing 1-bits in the binary representation of n. Note that tau(n) = v_2(n + 1) for odd n.

The stochastic model treats each iterate's parity as an independent Bernoulli(1/2) random variable. In the Terras formulation, the map T(n) = (3n + 1)/2 for odd n and T(n) = n/2 for even n combines the odd step and one mandatory even step into a single operation. Under this formulation, each operation changes log_2(n) by either +log_2(3/2) = +0.585 (odd) or -1 (even), each with probability 1/2, giving a net drift of -0.207 per operation.

In the uncompressed formulation C(n) used here, odd and even steps are counted separately. Since each odd step (3n + 1) is immediately followed by at least one even step (n/2), the fraction of odd steps among all steps is approximately 1/3 in the stochastic model. The net drift per individual step is then

Delta = (1/3) log_2(3) + (2/3)(-1) = 0.528 - 0.667 = -0.138 bits per step.

3. Validation: Per-Bit Stopping Time Increment

3.1 The Trailing-Bit Structure

A positive integer n with tau(n) = k trailing 1-bits satisfies n = 2^k * m - 1 for some integer m, or equivalently n ≡ 2^k - 1 (mod 2^k). The first k iterations of C on such a number are forced: each encounters an odd value, applies 3n + 1, then divides by 2 at least once. This contributes 2 forced steps per trailing bit (one odd, one even).

Each forced odd-then-even pair multiplies the value by approximately 3/2, increasing log_2(n) by 0.585 bits. To return this magnitude excess to equilibrium, the subsequent stochastic trajectory must shed 0.585 additional bits per trailing 1-bit. At the net drift rate of 0.138 bits per step, this requires 0.585 / 0.138 = 4.23 additional steps.

The predicted total increment per trailing 1-bit is therefore:

predicted = 2.00 (forced) + 4.23 (recovery) = 6.23 steps/bit.

3.2 Empirical Result

Computing mean stopping times grouped by trailing 1-bit count over n in [2, 10^7]:

The mean increment over tau = 2 through 5 (the stable regime with large sample sizes) is 6.17 steps/bit. The predicted value is 6.23 steps/bit. The discrepancy is 1.0%, a strong quantitative validation of the stochastic model's drift calculation.

3.3 Symmetry with 2-adic Valuation

The complementary analysis groups integers by v_2(n), the number of trailing 0-bits. Each trailing 0 reduces sigma(n) by one deterministic division, and the stochastic recovery cost is symmetric:

The decrement per trailing 0-bit (6.18 averaged over v_2 = 1 through 5) matches the increment per trailing 1-bit to within 0.2%, confirming the symmetry of the stochastic drift model.

3.4 Odd-Step Fraction

Among 100,000 random seeds in [2, 10^7], the mean fraction of odd steps is 0.3220, with standard deviation 0.0290. This is 3.4% below the theoretical 1/3, consistent with the known finite-size bias: smaller numbers (encountered as trajectories approach 1) have higher 2-adic valuation on average, increasing the fraction of even steps. This bias also explains why the observed 6.17 steps/bit is slightly below the predicted 6.23: the effective odd-step fraction in actual trajectories is lower than 1/3, producing a steeper net drift and fewer recovery steps.

4. Autocorrelation Structure

4.1 The Stopping-Time Landscape

The sequence {sigma(n)}_{n=1}^{N} is not a sequence of independent random variables. Adjacent integers can share trajectory fragments, and the identity

sigma(2m) = 1 + sigma(m) (*)

holds exactly for all m >= 1, creating a self-similar structure at all scales.

4.2 Autocorrelation Function

Define the autocorrelation at lag h as r(h) = Cov(sigma(n), sigma(n + h)) / Var(sigma(n)), computed over n in [2, N - h]. The empirical autocorrelation function over N = 10^7:

Two features are immediately apparent. First, r(1) = 0.402 is remarkably high: nearly half the variance in sigma(n) is shared with sigma(n + 1). Second, even at lag 1000, the autocorrelation remains at 0.067, far from zero.

4.3 Power-of-2 Lags and Self-Similarity

The identity (*) implies that sigma(n) and sigma(n + 2h) share more structure when h is a power of 2. Measuring autocorrelation exclusively at power-of-2 lags reveals a strikingly clean decay:

Each doubling of the lag reduces r by approximately 0.022. This near-linear relationship between r and k = log_2(h) implies logarithmic decay:

r(h) ~ 0.320 - 0.032 ln(h)

which can equivalently be expressed as a power law r(h) ~ 0.321 h^{-0.141}. This power-law fit achieves R^2 = 0.884 on the power-of-2 lags, vastly outperforming the exponential model r(h) ~ A exp(-h/tau), which yields R^2 = 0.315.

The extrapolated zero-crossing occurs at lag approximately 2^{14.5} ~ 22,000, suggesting that stopping times remain detectably correlated over a range of roughly 10^4 integers.

4.4 Even-Odd Lag Oscillation

The autocorrelation exhibits a pronounced even-odd oscillation: r(h) is systematically higher at even lags than at adjacent odd lags. For example, r(4) = 0.252 > r(3) = 0.212 and r(8) = 0.217 > r(7) = 0.152. This arises directly from (*): when n is even, sigma(n) = 1 + sigma(n/2), so pairs of even integers at lag 2h are structurally coupled through their half-values at lag h. Odd-lag correlations lack this recursive boost.

4.5 Structural Origin

The identity sigma(2m) = 1 + sigma(m) creates a fractal folding of the stopping-time sequence. Restricted to even indices, the sequence {sigma(2m)} is a vertically shifted copy of {sigma(m)}, compressed by a factor of 2 horizontally. This produces self-similar autocorrelation: the lag-2 autocorrelation among even seeds equals the lag-1 autocorrelation in the half-range population. We verified this identity empirically: both equal 0.3960 to all computed decimal places.

This recursive structure means the autocorrelation function is not determined by a single characteristic scale (as an exponential model would require), but inherits contributions from all dyadic scales simultaneously, producing the observed power-law decay.

5. Conditional Independence Tests

5.1 Variance Explained by Residue Classes

The stopping time sigma(n) depends on n mod 2^k through the first k forced parity steps. We measure the fraction of variance in sigma explained by conditioning on the residue class mod 2^k:

The variance explained grows at exactly 1.00% per bit, with no sign of saturation through k = 12. This perfect linearity means that each binary digit contributes equally to the systematic component of sigma(n), and that 88% of the variance at k = 12 is attributable to structure in the higher bits.

5.2 Deterministic Steps from Known Bits

To understand the per-bit contribution mechanistically, we measure how many trajectory steps are fully determined by the lowest k bits. Two seeds that agree in their lowest k bits but differ at bit k + 1 will follow identical Collatz trajectories until their parities first diverge. Over 5,000 random seeds in [10^6, 10^7]:

The relationship is almost perfectly linear: each known bit provides exactly 1.50 deterministic trajectory steps (R^2 = 0.999999, slope = 1.498). The ratio 3/2 arises because a known bit determines one parity test, and if the result is odd, the 3n + 1 step is followed by a guaranteed even step, adding 2 steps per odd encounter versus 1 step per even encounter, averaging 3/2 steps per bit.

5.3 Within-Class Autocorrelation

The central question of conditional independence: after grouping integers by their residue mod 2^k, is the within-class lag-1 autocorrelation zero? That is, for consecutive members of the class {n : n ≡ r (mod 2^k)}, namely r, r + 2^k, r + 2·2^k, ..., are their stopping times independent?

The within-class autocorrelation decreases with k but remains substantially positive at k = 6. Even after perfectly controlling for the lowest 6 bits, consecutive members of the same residue class have stopping times correlated at r = 0.106. This means that the binary suffix does not capture all sequential structure: the identity (*) and the correlations it induces operate at all scales, including scales far above 2^6.

The residual autocorrelation at lag 1000 within mod-64 classes drops to 0.050, compared to 0.067 in the unconditional sequence — a modest 25% reduction. The conditioning removes some short-range structure but leaves the long-range tail essentially intact.

6. Supplementary Results

6.1 Growth of Maximal Stopping Times

The ratio max_{n <= N} sigma(n) / (ln N)^2 provides an empirical test of the stochastic model's prediction for extremal behavior. The model of Kontorovich and Lagarias (2009) predicts that max sigma(n) grows as O((log N)^2), with the constant related to the extremal trajectory shape.

The ratio fluctuates around 2.6–2.7 for N >= 10^5, without clear convergence to a limit. The non-monotonicity (the ratio rises from 2.64 at N = 10^5 to 2.75 at N = 10^6 before returning to 2.64 at N = 10^7) reflects the sporadic nature of stopping-time records: new records can appear in clusters.

6.2 Record Stopping Times

There are 53 stopping-time records in [1, 10^7]. The record density (records per ln N) increases from 2.39 at N = 10^3 to 3.29 at N = 10^7, suggesting mildly superlogarithmic growth of the record count.

6.3 Distributional Shape

Stopping times are approximately log-normally distributed, with mean 155.27, standard deviation 61.76, and a slight negative skew of -0.371 in the log-transform. The Jarque-Bera statistic rules out exact log-normality (p ~ 0), indicating thin tails relative to the log-normal model (excess kurtosis = -0.188).

7. Discussion

7.1 The Stochastic Model Validated — and Bounded

The per-bit increment analysis (Section 3) provides the most precise test of the Terras-Lagarias stochastic model of which we are aware. The 1.0% agreement between the predicted 6.23 steps/bit and the observed 6.17 steps/bit demonstrates that the independent-parity assumption captures the average drift with high accuracy. The small residual (0.06 steps/bit) is accounted for by the finite-size bias in the odd-step fraction.

7.2 Long-Range Dependence

The autocorrelation analysis (Section 4) reveals structure that the independent-parity model cannot capture. The power-law decay r(h) ~ h^{-0.141} implies that the autocorrelation function is not integrable: the sum of r(h) over all h diverges. In time-series analysis, this is the hallmark of long-range dependence (also called long memory). The source is the exact arithmetic identity sigma(2m) = 1 + sigma(m), which recursively couples the stopping-time sequence at all dyadic scales.

The power-of-2 lag sequence provides the cleanest window into this structure. The autocorrelation decays by approximately 0.022 per doubling of the lag, corresponding to a loss of roughly 5.5% of the remaining correlation per doubling. This slow, self-similar decay is a signature of the binary tree structure underlying the Collatz iteration.

7.3 Conditional Independence and the Limits of Binary Conditioning

The conditional independence analysis (Section 5) shows that knowing the lowest k bits of n explains exactly k percent of the variance in sigma(n), but the within-class sequential dependence remains positive for all k tested. This means the autocorrelation structure arises not merely from shared low bits but from the global recursive structure of the Collatz map.

The relationship between known bits and deterministic steps (1.50 steps per bit) provides a mechanistic explanation for the 1% per-bit variance reduction: each bit pins down 1.50 trajectory steps, which at the overall trajectory length of ~155 steps represents approximately 1% of the total.

7.4 Open Questions

Several questions emerge from this analysis:

1. Does the autocorrelation exponent alpha = 0.141 have a closed-form expression? The proximity to log_2(3) - 1 = 0.585 or to 1/(2 log_2(3)) = 0.315 is not close enough to suggest an obvious formula, but the recursive structure of (*) should in principle determine alpha exactly.

2. Does the within-class autocorrelation converge to zero as k -> infinity, or does it approach a positive limit? Our data (k <= 6) cannot distinguish these possibilities. If the limit is positive, it would imply that no finite amount of binary conditioning can render stopping times independent — a fundamental obstruction to the independent-parity model.

3. The max sigma(n) / (ln N)^2 ratio has not clearly converged at N = 10^7. Data at N = 10^{12} or beyond, as computed by Oliveira e Silva and others, would be needed to establish the limiting constant.

4. The exact odd-step fraction of 0.322 (versus the theoretical 1/3) deserves a precise finite-size correction formula, potentially derivable from the known distribution of 2-adic valuations encountered along trajectories.

8. Conclusion

We have tested the standard stochastic model for Collatz stopping times against three quantitative predictions. The model passes the per-bit increment test with 1.0% accuracy, confirming that the independent-parity drift calculation correctly predicts average stopping time structure. However, the model fails to capture the strong sequential dependence between adjacent stopping times: the autocorrelation function follows a power law with exponent alpha ~ 0.14, indicating long-range memory in the stopping-time landscape. This dependence originates in the exact identity sigma(2m) = 1 + sigma(m), which the stochastic model cannot represent.

Conditioning on n mod 2^k removes variance at a rate of 1.00% per bit but does not eliminate the autocorrelation. The binary suffix captures the forced parity structure but not the global recursive coupling. These results delineate the precise boundary between what the stochastic model explains (average drift, per-bit increments) and what it cannot (sequential dependence, long-range correlation).

References

1. Terras, R. (1976). A stopping time problem on the positive integers. Acta Arithmetica, 30, 241–252.
2. Lagarias, J. C. (1985). The 3x + 1 problem and its generalizations. American Mathematical Monthly, 92(1), 3–23.
3. Kontorovich, A. V. and Lagarias, J. C. (2009). Stochastic models for the 3x + 1 and 5x + 1 problems. In The Ultimate Challenge: The 3x + 1 Problem, AMS, 131–188.
4. Tao, T. (2019). Almost all Collatz orbits attain almost bounded values. arXiv:1909.03562.
5. Lagarias, J. C. (2010). The 3x + 1 problem: An overview. In The Ultimate Challenge: The 3x + 1 Problem, AMS, 3–29.
6. Oliveira e Silva, T. (2010). Empirical verification of the 3x + 1 and related conjectures. In The Ultimate Challenge: The 3x + 1 Problem, AMS, 189–207.
7. Crandall, R. E. (1978). On the "3x + 1" problem. Mathematics of Computation, 32(144), 1281–1292.

Appendix A: Computational Details

All computations use Python 3 with only standard library functions. The random number generator is seeded with random.seed(42). Stopping times are computed by direct iteration of C(n) for all n in [1, 10^7], requiring approximately 45 seconds on a modern processor. Autocorrelation at lag h is computed as

r(h) = (1/M) sum_{n=2}^{N-h} (sigma(n) - mu)(sigma(n+h) - mu) / sigma^2

where mu and sigma^2 are the sample mean and variance over [2, N], and M = N - h - 1. For the conditional independence test, within-class autocorrelation at mod 2^k is the weighted average of lag-1 autocorrelations computed separately within each of the 2^k residue classes.