---
name: collatz-orbit-statistics
description: >
  Verifies the Collatz conjecture for all integers 1..10^6 and characterizes orbit
  statistics. For each n, computes stopping time (steps to reach 1) and max altitude
  (peak value in orbit) using a memoized Collatz iteration. Confirms all 1,000,000
  integers reach 1; measures mean stopping time growth as C × log₂(n) with C ≈ 6.59;
  identifies 44 delay-record integers; verifies known values: n=27 takes 111 steps,
  n=9663 reaches altitude 27,114,424, n=837799 takes 524 steps (longest in range).
  Zero pip installs, zero network, fully deterministic pure integer arithmetic.
  Triggers: Collatz conjecture, 3x+1 problem, orbit statistics, stopping time, delay
  records, hailstone sequence, Collatz verification, integer sequence benchmark.
allowed-tools: Bash(python3 *), Bash(mkdir *), Bash(cat *), Bash(cd *)
---

# Collatz Orbit Statistics

Verifies the Collatz conjecture for all positive integers up to 1,000,000 and measures
orbit statistics: stopping time (steps to reach 1), max altitude (peak orbit value),
delay records (numbers with stopping time exceeding all smaller numbers), and the empirical
constant C in the approximate growth law mean_stopping_time ≈ C × log₂(n).

The Collatz iteration: given n, apply n → n/2 if n is even, n → 3n+1 if n is odd, until
reaching 1. The conjecture states this always terminates; it has been verified computationally
up to approximately 2^68.

---

## Step 1: Setup Workspace

```bash
mkdir -p workspace && cd workspace
mkdir -p scripts output
```

Expected output:
```
(no terminal output — directories created silently)
```

---

## Step 2: Compute Collatz Orbits for n = 1..10^6

```bash
cd workspace
cat > scripts/collatz.py <<'PY'
#!/usr/bin/env python3
"""Collatz orbit statistics for n = 1..1,000,000.

Computes stopping time and max altitude for every integer in [1, 10^6]
using a memoized (cached) Collatz iteration, then reports:
  - Global summary statistics
  - Delay records (stopping-time record-holders)
  - Known-value cross-checks
  - Empirical C coefficient in mean ≈ C * log2(N)

Zero external dependencies. Pure integer arithmetic. Deterministic.
"""
import json
import math
import statistics

# ── Configurable parameter ────────────────────────────────────────────────────
N = 1_000_000
OUTPUT_FILE = "output/collatz_results.json"


def compute_collatz(n_max):
    """Compute stopping_time and max_altitude for all n in [1, n_max].

    Uses memoization: once stopping_time[x] is known, any orbit that
    reaches x can be resolved immediately without further iteration.

    Returns:
        stopping_time : list[int], length n_max+1, index 0 unused
        max_altitude  : list[int], length n_max+1, index 0 unused
    """
    stopping_time = [0] * (n_max + 1)
    max_altitude = [0] * (n_max + 1)

    # Base case: n=1 is already at 1
    stopping_time[1] = 0
    max_altitude[1] = 1

    for n in range(2, n_max + 1):
        # Walk the sequence until we hit 1 or a cached value
        path = []
        x = n
        while x != 1 and (x > n_max or stopping_time[x] == 0):
            path.append(x)
            x = x // 2 if x % 2 == 0 else 3 * x + 1

        # x is now a known anchor (either 1 or a cached value <= n_max)
        if x == 1:
            base_steps = 0
            base_max = 1
        else:
            base_steps = stopping_time[x]
            base_max = max_altitude[x]

        # Fill in the path from back to front
        cur_max = base_max
        steps_from_anchor = base_steps
        for val in reversed(path):
            steps_from_anchor += 1
            if val > cur_max:
                cur_max = val
            if val <= n_max:
                stopping_time[val] = steps_from_anchor
                max_altitude[val] = cur_max

    return stopping_time, max_altitude


def find_delay_records(stopping_time, n_max):
    """Return list of (n, steps) where stopping_time[n] > stopping_time[k] for all k < n."""
    records = []
    current_max = -1
    for n in range(1, n_max + 1):
        if stopping_time[n] > current_max:
            current_max = stopping_time[n]
            records.append({"n": n, "steps": stopping_time[n]})
    return records


def main():
    print(f"Computing Collatz orbits for n = 1..{N:,}...")
    stopping_time, max_altitude = compute_collatz(N)
    print("Computation complete.")

    # ── Verification: all n reach 1 ──────────────────────────────────────────
    all_reach_one = all(stopping_time[n] >= 0 for n in range(1, N + 1))
    # stopping_time[1]=0 and all others > 0 means orbit terminated
    all_positive = all(stopping_time[n] > 0 for n in range(2, N + 1))
    print(f"All n in [1, {N:,}] reach 1: {all_reach_one and all_positive}")

    # ── Global statistics ────────────────────────────────────────────────────
    stop_times = [stopping_time[n] for n in range(1, N + 1)]
    mean_st = statistics.mean(stop_times)
    stdev_st = statistics.stdev(stop_times)

    max_st = max(stop_times)
    max_st_n = stop_times.index(max_st) + 1  # offset by 1 (index 0 = n=1)
    max_alt = max(max_altitude[1:N + 1])
    max_alt_n = max_altitude.index(max_alt, 1)

    # ── C coefficient: mean / log2(N) ────────────────────────────────────────
    c_coeff = mean_st / math.log2(N)

    # ── Delay records ────────────────────────────────────────────────────────
    delay_records = find_delay_records(stopping_time, N)

    # ── Known-value cross-checks ─────────────────────────────────────────────
    known = {
        "n27_stopping_time":    stopping_time[27],
        "n27_max_altitude":     max_altitude[27],
        "n9663_stopping_time":  stopping_time[9663],
        "n9663_max_altitude":   max_altitude[9663],
        "n837799_stopping_time": stopping_time[837799],
    }

    # ── Assemble output ──────────────────────────────────────────────────────
    output = {
        "n_tested": N,
        "all_reach_one": all_reach_one and all_positive,
        "summary_stats": {
            "mean_stopping_time":   round(mean_st, 6),
            "stdev_stopping_time":  round(stdev_st, 6),
            "max_stopping_time":    max_st,
            "max_stopping_time_n":  max_st_n,
            "max_altitude":         max_alt,
            "max_altitude_n":       max_alt_n,
            "num_delay_records":    len(delay_records),
            "c_coefficient":        round(c_coeff, 6),
        },
        "known_values": known,
        "delay_records": delay_records,
    }

    with open(OUTPUT_FILE, "w") as fh:
        json.dump(output, fh, indent=2)

    # ── Print summary ─────────────────────────────────────────────────────────
    s = output["summary_stats"]
    print(f"n_tested            : {N:,}")
    print(f"all_reach_one       : {output['all_reach_one']}")
    print(f"mean_stopping_time  : {s['mean_stopping_time']}")
    print(f"stdev_stopping_time : {s['stdev_stopping_time']}")
    print(f"max_stopping_time   : {s['max_stopping_time']}  (n={s['max_stopping_time_n']})")
    print(f"max_altitude        : {s['max_altitude']}  (n={s['max_altitude_n']})")
    print(f"num_delay_records   : {s['num_delay_records']}")
    print(f"c_coefficient       : {s['c_coefficient']}  (mean / log2({N}))")
    print(f"n=27  steps={known['n27_stopping_time']}  max_alt={known['n27_max_altitude']}")
    print(f"n=9663 steps={known['n9663_stopping_time']}  max_alt={known['n9663_max_altitude']}")
    print(f"n=837799 steps={known['n837799_stopping_time']}")
    print(f"Results written to {OUTPUT_FILE}")


if __name__ == "__main__":
    main()
PY
python3 scripts/collatz.py
```

Expected output:
```
Computing Collatz orbits for n = 1..1,000,000...
Computation complete.
All n in [1, 1,000,000] reach 1: True
n_tested            : 1,000,000
all_reach_one       : True
mean_stopping_time  : 131.434424
stdev_stopping_time : 56.670087
max_stopping_time   : 524  (n=837799)
max_altitude        : 56991483520  (n=704511)
num_delay_records   : 44
c_coefficient       : 6.594284  (mean / log2(1000000))
n=27  steps=111  max_alt=9232
n=9663 steps=184  max_alt=27114424
n=837799 steps=524
Results written to output/collatz_results.json
```

---

## Step 3: Run Smoke Tests

```bash
cd workspace
python3 - <<'PY'
"""Smoke tests for Collatz orbit statistics."""
import json

results = json.load(open("output/collatz_results.json"))
s = results["summary_stats"]
kv = results["known_values"]

# ── Test 1: Exactly 1,000,000 starting values tested ──────────────────────────
assert results["n_tested"] == 1_000_000, \
    f"Expected n_tested=1000000, got {results['n_tested']}"
print("PASS  Test 1: exactly 1,000,000 starting values tested")

# ── Test 2: All reached 1 ─────────────────────────────────────────────────────
assert results["all_reach_one"] is True, \
    "FAIL: all_reach_one is not True"
print("PASS  Test 2: all 1,000,000 integers reach 1")

# ── Test 3: Mean stopping time is positive ────────────────────────────────────
assert s["mean_stopping_time"] > 0, \
    f"Mean stopping time should be positive, got {s['mean_stopping_time']}"
print(f"PASS  Test 3: mean stopping time > 0  (got {s['mean_stopping_time']:.4f})")

# ── Test 4: At least 10 delay records ─────────────────────────────────────────
assert s["num_delay_records"] >= 10, \
    f"Expected >= 10 delay records, got {s['num_delay_records']}"
print(f"PASS  Test 4: delay records list has >= 10 entries  (got {s['num_delay_records']})")

# ── Test 5: n=27 stopping time is exactly 111 ─────────────────────────────────
assert kv["n27_stopping_time"] == 111, \
    f"n=27 stopping time should be 111, got {kv['n27_stopping_time']}"
print("PASS  Test 5: n=27 stopping time is exactly 111")

# ── Test 6: n=837799 stopping time is exactly 524 ─────────────────────────────
assert kv["n837799_stopping_time"] == 524, \
    f"n=837799 stopping time should be 524, got {kv['n837799_stopping_time']}"
print("PASS  Test 6: n=837799 stopping time is exactly 524")

# ── Test 7: n=9663 max altitude is exactly 27,114,424 ─────────────────────────
assert kv["n9663_max_altitude"] == 27_114_424, \
    f"n=9663 max altitude should be 27114424, got {kv['n9663_max_altitude']}"
print("PASS  Test 7: n=9663 max altitude is exactly 27,114,424")

print()
print("smoke_tests_passed")
PY
```

Expected output:
```
PASS  Test 1: exactly 1,000,000 starting values tested
PASS  Test 2: all 1,000,000 integers reach 1
PASS  Test 3: mean stopping time > 0  (got 131.4344)
PASS  Test 4: delay records list has >= 10 entries  (got 44)
PASS  Test 5: n=27 stopping time is exactly 111
PASS  Test 6: n=837799 stopping time is exactly 524
PASS  Test 7: n=9663 max altitude is exactly 27,114,424

smoke_tests_passed
```

---

## Step 4: Print Delay Records Table

```bash
cd workspace
python3 - <<'PY'
"""Display delay records and C-coefficient growth table."""
import json
import math

results = json.load(open("output/collatz_results.json"))
records = results["delay_records"]

print("Delay records in [1, 1,000,000]")
print("(n whose stopping time exceeds all smaller n)")
print(f"{'n':>10}  {'steps':>6}")
print("-" * 20)
for r in records:
    print(f"{r['n']:>10}  {r['steps']:>6}")

print(f"\nTotal delay records: {len(records)}")
print()

# C coefficient growth table
s = results["summary_stats"]
print(f"Empirical C coefficient at N=1,000,000: {s['c_coefficient']}")
print("(mean_stopping_time / log2(N), theoretical asymptote ≈ 6.95)")
print(f"max_stopping_time = {s['max_stopping_time']}  (n={s['max_stopping_time_n']})")
print(f"max_altitude      = {s['max_altitude']}  (n={s['max_altitude_n']})")
PY
```

Expected output:
```
Delay records in [1, 1,000,000]
(n whose stopping time exceeds all smaller n)
         n   steps
--------------------
         1       0
         2       1
         3       7
         6       8
         7      16
         9      19
        18      20
        25      23
        27     111
        54     112
        73     115
        97     118
       129     121
       171     124
       231     127
       313     130
       327     143
       649     144
       703     170
       871     178
      1161     181
      2223     182
      2463     208
      2919     216
      3711     237
      6171     261
     10971     267
     13255     275
     17647     278
     23529     281
     26623     307
     34239     310
     35655     323
     52527     339
     77031     350
    106239     353
    142587     374
    156159     382
    216367     385
    230631     442
    410011     448
    511935     469
    626331     508
    837799     524

Total delay records: 44

Empirical C coefficient at N=1,000,000: 6.594284
(mean_stopping_time / log2(N), theoretical asymptote ≈ 6.95)
max_stopping_time = 524  (n=837799)
max_altitude      = 56991483520  (n=704511)
```

---

## Step 5: Verify Results

```bash
cd workspace
python3 - <<'PY'
import json

results = json.load(open("output/collatz_results.json"))
s = results["summary_stats"]
kv = results["known_values"]

# ── Core conjecture assertion ─────────────────────────────────────────────────
assert results["all_reach_one"] is True, \
    "Conjecture verification FAILED: not all integers in [1, 10^6] reach 1"

# ── Known-answer assertions ───────────────────────────────────────────────────
assert kv["n27_stopping_time"] == 111, \
    f"n=27 should take 111 steps, got {kv['n27_stopping_time']}"

assert kv["n9663_max_altitude"] == 27_114_424, \
    f"n=9663 max altitude should be 27114424, got {kv['n9663_max_altitude']}"

assert s["max_stopping_time_n"] == 837_799, \
    f"Longest stopping time should be n=837799, got n={s['max_stopping_time_n']}"

assert s["max_stopping_time"] == 524, \
    f"Longest stopping time should be 524 steps, got {s['max_stopping_time']}"

# ── Statistical plausibility assertions ──────────────────────────────────────
assert 120 < s["mean_stopping_time"] < 145, \
    f"Mean stopping time out of expected range: {s['mean_stopping_time']}"

assert 5.5 < s["c_coefficient"] < 7.5, \
    f"C coefficient out of expected range: {s['c_coefficient']}"

assert s["num_delay_records"] == 44, \
    f"Expected 44 delay records, got {s['num_delay_records']}"

assert s["max_altitude_n"] == 704_511, \
    f"Highest altitude should be n=704511, got n={s['max_altitude_n']}"

print("Collatz conjecture verified for all n in [1, 1,000,000].")
print(f"  Longest orbit : n=837799  (524 steps)")
print(f"  Highest peak  : n=704511  (altitude 56,991,483,520)")
print(f"  Delay records : 44 integers in [1, 10^6]")
print(f"  C coefficient : {s['c_coefficient']}  (theoretical asymptote ≈ 6.95)")
print()
print("collatz_orbit_statistics_verified")
PY
```

Expected output:
```
Collatz conjecture verified for all n in [1, 1,000,000].
  Longest orbit : n=837799  (524 steps)
  Highest peak  : n=704511  (altitude 56,991,483,520)
  Delay records : 44 integers in [1, 10^6]
  C coefficient : 6.594284  (theoretical asymptote ≈ 6.95)

collatz_orbit_statistics_verified
```

---

## Notes

### What This Measures

The **stopping time** T(n) is the number of steps for the Collatz iteration to reach 1.
T(1)=0, T(2)=1, T(3)=7, ..., T(27)=111, T(837799)=524.

The **max altitude** A(n) is the maximum value ever reached in the orbit of n,
including n itself. A(27)=9,232 despite T(27)=111.

**Delay records** are integers n satisfying T(n) > T(k) for all k < n. They mark the
frontier of orbit complexity as n grows. In [1, 10^6] there are exactly 44 such records,
starting at {1, 2, 3, 6, 7, 9, ...} and ending at 837799.

The **C coefficient** measures how mean stopping time grows with n: empirically
C(N) = mean(T(1..N)) / log₂(N) ≈ 6.59 at N=10^6, approaching the theoretical
asymptote C∞ ≈ 6.95 as N → ∞.

### Memoization

The algorithm walks each orbit forward until reaching either 1 (the base case) or
a previously computed value. It then fills in the entire un-cached prefix in one pass.
This reduces the total number of Collatz steps computed from O(N × mean_steps) ≈ 130 × 10^6
to roughly O(N × log N) in practice, cutting runtime to under 30 seconds in Python.

### Limitations

1. **10^6 is a tiny verification window.** The conjecture has been verified up to approximately
   2^68 ≈ 2.95 × 10^20 (Roosendaal, ongoing). This skill covers only a minuscule fraction
   of that frontier.

2. **Orbit statistics do not prove the conjecture.** Measuring mean stopping time and delay
   records is descriptive; no finite verification can exclude counterexamples at arbitrarily
   large n.

3. **The C coefficient has not converged.** C(10^6) ≈ 6.59 versus the theoretical 6.95 —
   a ~5% underestimate. Extending to N = 10^9 would narrow this gap.

4. **Max altitude grows faster than stopping time.** A(704511) ≈ 5.7 × 10^10 despite
   T(704511) being only moderate; altitude extremes are harder to characterize statistically
   than stopping times.

5. **No log-normal fit computed.** The stopping-time distribution is approximately log-normal
   but this skill only reports mean and standard deviation; a formal Kolmogorov-Smirnov test
   against the log-normal CDF is left for future work.