Introduction

The Collatz map sends a positive integer \(n\) to \(n/2\) if it is even and to \(3n + 1\) if it is odd; the (still unproven) conjecture is that every starting value eventually reaches \(1\). All four figures below are different renderings of the same object — either the reverse (predecessor) tree of the map, or a single long orbit — and everything here is pure base R, so the only package you need is rmarkdown to knit.

Two figures are deterministic (fixed geometry / a fixed orbit) and reproduce the originals essentially exactly:

  • Figure 2 — a direct port of the Processing “Collatz Tree” sketch.
  • Figure 4 — the polar Proportional Power Ratio web of one long orbit.

The other two come from graph layouts. The originals used graphviz neato and a D3 force simulation; here I use a self-contained radial-tree layout and a vectorised Fruchterman–Reingold layout. The structure and colouring match, but the exact placement is layout- and seed-dependent, so they are faithful stylistic reproductions rather than pixel-identical copies. Change the seed argument to reshuffle them.

## ---- Shared Collatz machinery -------------------------------------------

# Forward "3n+1" orbit until it reaches 1
collatz_orbit <- function(n, max_steps = 2000) {
  n <- as.double(n); seq <- n
  while (n != 1 && length(seq) < max_steps) {
    n <- if (n %% 2 == 0) n / 2 else 3 * n + 1
    seq <- c(seq, n)
  }
  seq
}

# Proportional Power Ratio of base b: maps x into [0,1) by its position
# between consecutive powers of b.  P(x) = (x - b^k)/(b^(k+1) - b^k)
ppr <- function(x, b = 6) {
  k  <- floor(log(x, base = b))
  bk <- b^k
  (x - bk) / (b * bk - bk)
}

# Reverse-Collatz predecessor tree (shortcut version, matches the Raku edges):
#   n -> 2n         (always)
#   n -> (2n-1)/3   (only when n %% 3 == 2)
build_reverse_tree <- function(root = 64, levels = 18, max_nodes = 2200) {
  val <- root; parent <- 0L; depth <- 0L
  frontier <- 1L
  for (lv in 1:levels) {
    newf <- integer(0)
    for (i in frontier) {
      n <- val[i]; kids <- 2 * n
      if (n %% 3 == 2) kids <- c(kids, (2 * n - 1) / 3)
      for (k in kids) {
        val <- c(val, k); parent <- c(parent, i); depth <- c(depth, lv)
        newf <- c(newf, length(val))
        if (length(val) >= max_nodes) break
      }
      if (length(val) >= max_nodes) break
    }
    frontier <- newf
    if (length(val) >= max_nodes || length(frontier) == 0) break
  }
  list(value = val, parent = parent, depth = depth)
}

Figure 1 — The green/red “neuron” (force-directed tree)

A large predecessor tree laid out with a Fruchterman–Reingold force simulation. Edges are drawn as quadratic Béziers, which produces the organic loops of the original; the longest root-to-leaf paths are highlighted in red so a few backbone strands read through the tangle.

# Vectorised Fruchterman-Reingold force layout
fr_layout <- function(from, to, n, iterations = 240, seed = 7, cool = 0.985) {
  set.seed(seed)
  pos  <- matrix(rnorm(2 * n, sd = 0.5), n, 2)
  k    <- sqrt(1 / n) * 2.2                     # ideal edge length
  temp <- 0.6
  for (it in 1:iterations) {
    dx <- outer(pos[,1], pos[,1], "-"); dy <- outer(pos[,2], pos[,2], "-")
    dist <- sqrt(dx*dx + dy*dy); dist[dist < 1e-3] <- 1e-3
    rep <- (k*k) / dist                         # repulsion
    fx <- rowSums(dx / dist * rep); fy <- rowSums(dy / dist * rep)
    ex <- pos[from,1] - pos[to,1]; ey <- pos[from,2] - pos[to,2]
    ed <- sqrt(ex*ex + ey*ey); ed[ed < 1e-3] <- 1e-3
    att <- (ed*ed) / k                          # attraction along edges
    ax <- ex / ed * att; ay <- ey / ed * att
    fx <- fx - tapply2(ax, from, n) + tapply2(ax, to, n)
    fy <- fy - tapply2(ay, from, n) + tapply2(ay, to, n)
    dl <- sqrt(fx*fx + fy*fy); dl[dl < 1e-6] <- 1e-6
    pos[,1] <- pos[,1] + (fx / dl) * pmin(dl, temp)
    pos[,2] <- pos[,2] + (fy / dl) * pmin(dl, temp)
    temp <- temp * cool
  }
  pos
}
tapply2 <- function(vals, idx, n) {             # grouped sum -> length-n vector
  out <- numeric(n); tab <- tapply(vals, idx, sum)
  out[as.integer(names(tab))] <- tab; out
}

# node indices whose edge-to-parent lies on one of the K longest paths
top_paths <- function(tree, K = 6) {
  n <- length(tree$value); pr <- tree$parent; d <- tree$depth
  children <- split(seq_len(n), pr)
  leaves <- setdiff(seq_len(n), as.integer(names(children)))
  leaves <- leaves[order(d[leaves], decreasing = TRUE)][seq_len(min(K, length(leaves)))]
  edges <- integer(0)
  for (lf in leaves) { i <- lf; while (pr[i] > 0) { edges <- c(edges, i); i <- pr[i] } }
  unique(edges)
}

# draw many quadratic-bezier curved edges at once
curved_segments <- function(x0, y0, x1, y1, col, lwd, bend = 0.28, npts = 16) {
  mx <- (x0 + x1)/2; my <- (y0 + y1)/2
  dx <- x1 - x0; dy <- y1 - y0
  cx <- mx - dy * bend; cy <- my + dx * bend    # control point (perp. offset)
  t  <- seq(0, 1, length.out = npts)
  px <- outer(x0, (1-t)^2) + outer(cx, 2*(1-t)*t) + outer(x1, t^2)
  py <- outer(y0, (1-t)^2) + outer(cy, 2*(1-t)*t) + outer(y1, t^2)
  for (j in 1:(npts-1)) segments(px[,j], py[,j], px[,j+1], py[,j+1], col = col, lwd = lwd)
}

plot_image1 <- function(root = 8, levels = 22, max_nodes = 1600,
                        iterations = 170, seed = 11, bend = 0.30) {
  tr <- build_reverse_tree(root, levels, max_nodes)
  n  <- length(tr$value); pr <- tr$parent; e <- which(pr > 0)
  pos <- fr_layout(from = e, to = pr[e], n = n, iterations = iterations, seed = seed)
  is_red <- e %in% top_paths(tr, K = 9)
  op <- par(bg = "black", mar = c(0,0,0,0)); on.exit(par(op))
  plot(pos, type = "n", asp = 1, axes = FALSE, xlab = "", ylab = "")
  gr <- e[!is_red]; rd <- e[is_red]
  curved_segments(pos[gr,1], pos[gr,2], pos[pr[gr],1], pos[pr[gr],2],
                  col = "#2FCF5A", lwd = 1.3, bend = bend)
  curved_segments(pos[rd,1], pos[rd,2], pos[pr[rd],1], pos[pr[rd],2],
                  col = "#E23A1E", lwd = 1.5, bend = bend)
  points(pos, pch = 16, col = "#1EA84B", cex = 0.42)
}
# ~1600 nodes + 170 FR iterations takes ~20 s; drop max_nodes/iterations to speed up
plot_image1(root = 8, levels = 22, max_nodes = 1600, iterations = 170, seed = 11)

Figure 2 — The rainbow fan (Processing “Collatz Tree” port)

A faithful port of the Processing sketch. Starting from the root num = 1 at angle \(\pi/4\), every node branches to \(2n\) (even branch, angle \(+a_e\)) and, when \((n-1)\bmod 3 = 0\), to \((n-1)/3\) (odd branch, angle \(+a_o\)). All branch lengths are equal, colours cycle through the original colray[] palette by depth, and leaves are drawn as green dots. This one is fully deterministic.

build_collatz_tree <- function(max_depth = 24, branchlen = 26,
                               ae = 4 * pi / 90, ao = -8 * pi / 90) {
  num <- 1; x <- 0; y <- 0; ang <- pi / 4
  seg_x1 <- seg_y1 <- seg_x2 <- seg_y2 <- seg_d <- numeric(0)
  for (depth in 0:(max_depth - 1)) {
    ea <- ang + ae
    ex <- x + branchlen * cos(ea); ey <- y + branchlen * sin(ea)
    seg_x1 <- c(seg_x1, x); seg_y1 <- c(seg_y1, y)
    seg_x2 <- c(seg_x2, ex); seg_y2 <- c(seg_y2, ey)
    seg_d  <- c(seg_d, rep(depth, length(x)))
    n2 <- num * 2; x2 <- ex; y2 <- ey; a2 <- ea
    keep <- which(num > 1 & (num - 1) %% 3 == 0)
    if (length(keep)) {
      oa <- ang[keep] + ao
      ox <- x[keep] + branchlen * cos(oa); oy <- y[keep] + branchlen * sin(oa)
      seg_x1 <- c(seg_x1, x[keep]); seg_y1 <- c(seg_y1, y[keep])
      seg_x2 <- c(seg_x2, ox); seg_y2 <- c(seg_y2, oy)
      seg_d  <- c(seg_d, rep(depth, length(keep)))
      n2 <- c(n2, (num[keep]-1)/3); x2 <- c(x2, ox); y2 <- c(y2, oy); a2 <- c(a2, oa)
    }
    num <- n2; x <- x2; y <- y2; ang <- a2
  }
  list(seg = data.frame(x1=seg_x1, y1=seg_y1, x2=seg_x2, y2=seg_y2, depth=seg_d),
       leaf = data.frame(x = x, y = y))
}

# rainbow palette taken from the Processing colray[] (0xffRRGGBB -> #RRGGBB)
collatz_palette <- local({
  argb <- c(0x0000ff,0x2000ff,0x4000ff,0x6000ff,0x8000ff,0xA000ff,0xC000ff,0xE000ff,
            0xff00E0,0xff00C0,0xff00A0,0xff0080,0xff0060,0xff0040,0xff6020,0xff8000,
            0xffA000,0xffC000,0xffE000,0xffff00,
            0xffE000,0xffC000,0xffA000,0xff8000,0xff6020,0xff0040,0xff0060,0xff0080,
            0xff00A0,0xff00C0,0xff00E0,0xE000ff,0xC000ff,0xA000ff,0x8000ff,0x6000ff,
            0x4000ff,0x2000ff)
  sprintf("#%06X", argb)
})

plot_image2 <- function(max_depth = 24) {
  tr <- build_collatz_tree(max_depth); s <- tr$seg
  op <- par(bg = "black", mar = c(0,0,0,0)); on.exit(par(op))
  plot(NA, xlim = range(c(s$x1, s$x2)), ylim = range(c(s$y1, s$y2)),
       asp = 1, axes = FALSE, xlab = "", ylab = "")
  segments(s$x1, s$y1, s$x2, s$y2,
           col = collatz_palette[(s$depth %% length(collatz_palette)) + 1], lwd = 2)
  points(tr$leaf$x, tr$leaf$y, pch = 16, col = "#00FF00", cex = 0.7)
}
plot_image2(max_depth = 24)

Figure 3 — The purple/orange neato-style tree

The same predecessor tree, but laid out with the Fruchterman–Reingold force simulation (reusing fr_layout from Figure 1). This is what gives the organic graphviz-neato look of the original: major subtrees physically repel each other, so the tree splays into distinct lobes separated by white gaps, with leaf-fans spreading at the rim. Edges are straight purple lines; nodes are coloured by depth from the root, following the Raku blend — dark-red → plum over the inner levels, then plum → orange out to the leaves.

The exact lobe arrangement is seed-dependent (a force simulation has no canonical solution); change seed to rotate/reshape the lobes.

plot_image3 <- function(root = 64, levels = 20, max_nodes = 1928,
                        iterations = 300, seed = 3) {
  tr  <- build_reverse_tree(root, levels, max_nodes)
  n   <- length(tr$value); pr <- tr$parent; d <- tr$depth
  e   <- which(pr > 0)
  pos <- fr_layout(from = e, to = pr[e], n = n, iterations = iterations, seed = seed)
  # depth palette, matching the Raku two-segment blend darkred -> plum(16) -> orange.
  # (Most nodes sit at large depth, so the plum->orange leg fills the visible rim.)
  nc  <- max(d) + 1
  pal <- c(colorRampPalette(c("#8B0000", "#DDA0DD"))(min(16, nc)),
           if (nc > 16) colorRampPalette(c("#DDA0DD", "#FFA500"))(nc - 16))
  vcol <- pal[d + 1]
  csz  <- 0.9 - 0.42 * (d / max(d))    # larger dark-red hubs, small orange leaves
  op <- par(bg = "white", mar = c(0,0,0,0)); on.exit(par(op))
  plot(pos, type = "n", asp = 1, axes = FALSE, xlab = "", ylab = "")
  segments(pos[e,1], pos[e,2], pos[pr[e],1], pos[pr[e],2], col = "#8E2C8E", lwd = 1.0)
  points(pos, pch = 16, col = vcol, cex = csz)
}
# ~1900 nodes x 300 FR iterations (~40 s). Lower max_nodes/iterations for a draft;
# change `seed` to rotate/reshape the lobes.
plot_image3(root = 64, levels = 20, max_nodes = 1928, iterations = 300, seed = 3)

Figure 4 — The polar Proportional-Power-Ratio web

A single long orbit, plotted in polar coordinates with \(\theta = \text{index}\) (in radians) and \(r = P(\text{value})\) in base 6. Because each step advances the angle by one radian, the orbit winds many times around the circle; connecting consecutive points draws the dense star/web, and the near-resonance of the index against \(2\pi\) produces the characteristic spokes. The default seed n = 670617279 has a 949-step orbit. This figure is exact.

plot_image4 <- function(n = 670617279, b = 6) {
  orb <- collatz_orbit(n)
  theta <- seq_along(orb) - 1
  r <- ppr(orb, b); X <- r * cos(theta); Y <- r * sin(theta)
  op <- par(bg = "white", mar = c(0,0,0,0)); on.exit(par(op))
  plot(NA, xlim = c(-1.15,1.15), ylim = c(-1.15,1.15),
       asp = 1, axes = FALSE, xlab = "", ylab = "")
  ang <- seq(0, 2*pi, length.out = 400)
  for (rr in c(0.6,0.8,1.0)) lines(rr*cos(ang), rr*sin(ang), col = "grey75", lwd = 0.8)
  lines(cos(ang), sin(ang), col = "grey30", lwd = 1)
  for (a in seq(0, 2*pi, by = pi/4)) segments(0,0, cos(a), sin(a), col = "grey85", lwd = 0.6)
  for (a in seq(0, 315, by = 45)) {
    ar <- a*pi/180; text(1.09*cos(ar), 1.09*sin(ar), paste0(a,"\u00B0"), col="grey40", cex=0.8)
  }
  for (rr in c(0.6,0.8,1.0)) text(rr*cos(pi/12), rr*sin(pi/12), format(rr), col="grey40", cex=0.7)
  lines(X, Y, col = "#1f77b4", lwd = 0.9)
}
plot_image4(n = 670617279, b = 6)

Notes & knobs

  • Deterministic vs. layout-dependent. Figures 2 and 4 are fixed geometry and reproduce the sources exactly. Figures 1 and 3 depend on a force / radial layout — same tree, same colouring, but placement shifts with seed, max_nodes, and iterations.
  • Figure 1 runtime. Fruchterman–Reingold here is \(O(n^2)\) per iteration, so ~1600 nodes × 170 iterations is ~20 s. Lower max_nodes or iterations for a quick draft; raise them for a denser blob.
  • Try other seeds/orbits. Interesting Figure-4 seeds include n = 27, n = 97, and n = 703. For Figure 2, sweep the branch angles ae/ao inside build_collatz_tree() to reshape the fan.
  • If you prefer igraph: Figures 1 and 3 map cleanly onto graph_from_edgelist() + layout_with_fr() / layout_as_tree(circular=TRUE) if you’d rather lean on a library — the base-R versions above just keep the file dependency-free. ```