the r in recursion stands for recursion

problem

define $s: \mathbb{N}_0^4 \mapsto \mathbb{N}_0^4; s(a, b, c, d) \mapsto \left(|a - b|, |b - c|, |c - d|, |d - a|\right)$, and let $f: \mathbb{N}_0^4 \mapsto \mathbb{N}_0$ be the number of times needed to recursively apply $s$ to $(a, b, c, d)$ to get $(0, 0, 0, 0)$; in other words, $s^{f(a, b, c, d)}(a, b, c, d) = (0, 0, 0, 0)$.

for example, $f(0, 0, 0, 0)$ is trivially $0$, while $f(6, 4, 2, 4) = 2$ since $s(6, 4, 2, 4) = (2, 2, 2, 2)$ and $s(2, 2, 2, 2) = (0, 0, 0, 0)$.

maximize $f$ over $0 \leq a, b, c, d \leq 10^7$, and state where such maxima occur.

source: jane street puzzles, jan 2023

---

solution

the first thing to observe is that $s$ has some symmetry: if $s(a, b, c, d) = (w, x, y, z)$ then $s(b, c, d, a) = (x, y, z, w)$. this implies that $f$ is cyclic: any “shifts” of $(a, b, c, d)$ will not change the value of $f$. to put this into a picture, imagine putting $a, b, c, d$ down on a circle, in that order, so you could do something like $a$ in the 12 o’clock position, $b$ at 3 o’clock, $c$ at 6, $d$ at 9. then, starting at any of the four points on this circle and going either clockwise or counterclockwise, $f$ of that permutation of $(a, b, c, d)$ you read off will equal $f(a, b, c, d)$. for example, if you started from $c$ and went counterclockwise you’d get $f(c, b, a, d) = f(a, b, c, d)$.

now we can simplify $(a, b, c, d)$. without loss of generality assume $d$ is the (distinctly) smallest; then since $s(a, b, c, d) = s(a - d, b - d, c - d, 0)$ we have $f(a, b, c, d) = f(a - d, b - d, c - d, 0)$. since $f(c, b, a, d) = f(a, b, c, d)$ we may also assume without loss of generality that $a > c$. (the cases of if one of $a, b, c = d$, or $a, b = c$, is just computation and has $f(a, b, c, d)$ at most like 10 something.)

so, we examine $f(a, b, c, 0)$ with $a > c$; we claim that $a > b > c$ maximizes $f$. let’s try the other two cases, then: if $b > a > c$ then $s(a, b, c, 0) = (b - a, b - c, c, a)$, and noticing $b - c = (b - a) + (a - c)$ we have $s(b - a, b - c, c, a) = \left(a - c, |b - 2c|, a - c, |b - 2a|\right)$. but now note that $s(x, y, x, z) = \left(|x - y|, |x - y|, |x - z|, |x - z|\right)$ and two more applications of $s$ will reduce it down to $(0, 0, 0, 0)$, so $f(a, b, c, 0)$ with $b > a > c$ is at most 4; similarly $f(a, b, c, d)$ with $a > c > b$ is also at most 4.

we now change the domain of $s$ and $f$; instead of $\mathbb{N}_0^4$ it’ll be $\mathbb{R}_{\geq 0}^4$, over the positive reals. this is so we can further take a degree of freedom away from $f(a, b, c, 0)$: dividing through by $a$ leaves $f(1, p, q, 0)$ for rationals $1 > p > q$. note that $s$ is continuous, so $f$ should be “semi-continuous”: $f$ returns integers and so is clearly not continuous, but the point is that small changes in $p, q$ should only result in $f$ going up or down by 1, the minimum unit of change. in other words, for small rational $\epsilon$, $f(1, p \pm \epsilon, q \pm \epsilon, 0)$ should differ from $f(1, p, q, 0)$ by no more than 1 in either direction.

it remains to see if there exists reals $p, q$ such that $s(1, p, q, 0) = (1, p, q, 0)$ (or some cyclic equivalent of that); this is called a fixed point. this is because then $f(1, p, q, 0) = \infty$, so all rationals $p’, q’$ close to $p, q$ will have very high $f$ values. the first step is to simplify $s(1, p, q, 0) = (1 - p, p - q, q, 1)$ into the form $(1, x, y, 0)$. to get the appropriate cyclic equivalent of $(1 - p, p - q, q, 1)$, we know the highest of these four must be first, so $1$ must be first. placing the four numbers $1 - p, p - q, q, 1$ onto the circle and starting off with 1, we have two cases (depending on if we go clockwise or counterclockwise): $1 > 1 - p > p - q > q$ or $1 > q > p - q > 1 - p$, which produce simplified pairs of $\left(1, \frac{1 - p - q}{1 - q}, \frac{p - 2q}{1 - q}, 0\right)$ and $\left(1, \frac{p + q - 1}{p}, \frac{2p - q - 1}{p}, 0\right)$ respectively.

the first case gives us $p(1 - q) = 1 - p - q$ and $q(1 - q) = p - 2q \iff p = q(3 - q)$; substituting back into the first equation gives us the cubic $q^3 - 5q^2 + 7q - 1$ which does have an appropriate root $q \in (0, 1)$. however, the second case gives us $p^2 = p + q - 1 \iff q = p^2 - p + 1$ and $pq = 2p - q - 1 \iff p^3 - 2p + 2$ which does not have an appropriate root $p \in (0, 1)$. so, let $(p, q)$ be the solution we get from our first case, $q^3 - 5q^2 + 7q - 1 = 0$.

it remains to find rational approximations to $p \approx \frac{b}{a}$ and $q \approx \frac{c}{a}$, getting solutions $(a, b, c, 0)$ from there, then calculating $f(a, b, c, 0)$ over all approximations to find the maximum. we can save some time by iterating over $c$, the smallest, instead of $a$, since then $a \approx \frac{c}{q}, b \approx ap = \frac{cp}{q}$. this means that we only need to iterate over $1 \leq c \leq \left\lceil 10^7q \right\rceil$ since $a \leq 10^7$. note that our code checks both the floor and ceiling of $\frac{c}{q}$ and $\frac{cp}{q}$ since it’s not clear which approximation will give the better results.

def s(a, b, c, d):
    return (abs(a - b), abs(b - c), abs(c - d), abs(d - a))

def f(a, b, c, d):
    if (a, b, c, d) == (0, 0, 0, 0): 
        return 1
    else:
        (a, b, c, d) = s(a, b, c, d)
        return 1 + f(a, b, c, d)

# plug into wolframalpha to actually get the root q
q = (5 - 4 / ((19 - 3 * (33 ** (1/2))) ** (1/3)) - ((19 - 3 * (33 ** (1/2))) ** (1/3))) / 3
p = q * (3 - q)
limit = 10 ** 7 + 1

max_f = 0
result = (limit, limit, limit, limit)

for c in range(math.ceil((10 ** 7) * q)): 
    a = math.floor(c / q)
    b = math.floor(a * p)
    # f1, f2, f3, f4 are all needed to check both floor and ceiling of our approximations
    f1 = f(a, b, c, 0)
    if f1 > max_f:
        max_f = f1
        result = (a, b, c, 0)
    f2 = f(a, b + 1, c, 0)
    if f2 > max_f:
        max_f = f2
        result = (a, b + 1, c, 0)
    f3 = f(a + 1, b, c, 0)
    if f3 > max_f:
        max_f = f3
        result = (a + 1, b, c, 0)
    f4 = f(a + 1, b + 1, c, 0)
    if f4 > max_f:
        max_f = f4
        result = (a + 1, b + 1, c, 0)
        
print(max_f)
print(result)

our maximum is then $f(8646064, 3945294, 1389537, 0) = 43.$ we are done.