Consistent Collatz

A consistent Collatz sequence is a sequence of integers which can be defined by a recurrence relation of the following form:

$mx_{n+1}=rx_{n}+J[x_{n}\ {\textrm {mod}}\ m]$

where the sequence itself is $x_{0},x_{1},\ldots$, $m$ and $r$ are positive integers, and $J$ is a map from integers modulo $m$ to (possibly negative) integers (with the values chosen to ensure that all $x_{n}$ are integers). The sequence is entirely defined by $m$, $r$, $J$, and the value of its first element $x_{0}$.

Many small cryptids work by calculating the elements of a consistent Collatz sequence modulo $m$, and deciding whether or not to halt based on the remainders. For example, Hydra calculates the consistent Collatz sequence with $m=2$, $r=3$, $J=\{0\mapsto 0,1\mapsto -1\}$, and $x_{0}=3$, halting only if the sequence has contained more than twice as many even elements as odd elements. $m=2,r=3$ seems to be particularly common among small cryptids, due to creating the simplest consistent Collatz sequences that have nontrivial behaviour; however, other values have been observed, such as Bigfoot's $m=81,r=256$.

Efficiently calculating consistent Collatz sequences

For any given Collatz sequence, it is possible to calculate its sequence of remainders $x_{n}\ {\textrm {mod}}\ m$ in amortized quasilinear time. This can be accomplished via the use of two helper sequences:

$w_{n}=x_{n}\ {\textrm {mod}}\ m^{f(n)}$, except $w_{0}=x_{0}$
$j_{n}=m^{f(n)}x_{n}-r^{f(n)}x_{n-f(n)}$, with $j_{0}$ undefined (the algorithm never uses it)

where $f(n)$ is the largest power of 2 that divides into $n$.

Time complexity

The algorithm works by calculating $j_{n}$ then $w_{n}$ for each $n$ in turn, with the only operations used being additions, subtractions, and multiplications of numbers whose number of digits is proportional to $f(n)$; division and modulo by powers of $m$; and calculation of $r^{t}$ with $t\leq f(n)$. Because the values of $r^{t}$ can be memoized, and it is possible to trivialise the division and modulus operations by storing the numbers in base $m$, this means that the only slow operations are additions and subtractions taking $O(f(n))$ time, and multiplications taking $O(f(n)\log f(n)\log \log f(n))$ time (and there are $O(\log f(n))$ such operations performed), so the calculation of each $w_{n},j_{n}$ pair takes time quasilinear in $f(n)$. $\Sigma _{t=1}^{n}{f(t)}=O(n\log n)$: as such, the calculation of the entire sequences $w_{0},\ldots ,w_{n}$ and $j_{0},\ldots ,j_{n}$ takes time quasilinear in $n\log n$, thus quasilinear in $n$.

Details of the algorithm

The algorithm itself is:

$j_{n}=m^{f(n)-1}J[w_{n-1}\ {\textrm {mod}}\ m]+\Sigma _{e=1}^{\log _{2}f(n)}\left(m^{f(n)-2^{e}}r^{2^{e-1}}j_{n-2^{e-1}}\right)$
$w'_{n}=j_{n}+r^{f(n)}w_{n-f(n)}\mod m^{2f(n)}$ (i.e. all the calculations are done modulo $m^{2f(n)}$, saving time in cases where $w_{n-f(n)}$ happens to be much larger than $m^{2f(n)}$)
$w_{n}=w'_{n}\div m^{f(n)}$ (which should always be an integer).

(If $n$ is odd, then $f(n)=1$, $\log _{2}f(n)=0$, and the sum in the calculation of $j_{n}$ is a degenerate sum with no elements.)

Sketch proof of correctness

The proof that the algorithm is correct starts with the definition $j_{n}=m^{f(n)}x_{n}-r^{f(n)}x_{n-f(n)}$, adds a degenerate sum (that sums to zero) to produce the following expression:

$j_{n}=m^{f(n)}x_{n}+\Sigma _{e=1}^{\log _{2}f(n)}\left(-m^{f(n)-2^{e-1}}r^{2^{e-1}}x_{n-2^{e-1}}+m^{f(n)-2^{e-1}}r^{2^{e-1}}x_{n-2^{e-1}}\right)-r^{f(n)}x_{n-f(n)}$ (with each element of the sum being of the form $-q+q$ and thus 0)

and then rebrackets by grouping the term before the sum with the first half of the first element of the sum, the second half of the first element of the sum with the first half of the second element of the sum, etc., and finally the second half of the last element of the sum with the term after the sum. From there, the proof is mostly just a matter of expanding definitions.

Possible tricks to optimize the implementation

Once the algorithm produces a value of $j_{n}$ or $w_{n}$, it happens that it will never again read a value $j_{n'}$ or $w_{n'}$ where $f(n)=f(n')$. As such, it is possible to save memory by storing only one value of $j_{n}$ and one value of $w_{n}$ for each $f(n)$.

It is probably faster to switch to an alternative implementation (e.g. simple repeated multiplication) when the $f(n)$ values are small, because at that point the numbers are small enough to fit into a machine register; this does not improve the asymptotic behaviour of the implementation but is likely to speed it up by a reasonably high constant factor.

Every $w'_{n}$ is necessarily a multiple of $m^{f(n)}$. It seems like that might provide some sort of shortcut to calculate it faster, although the details are currently unclear.

Proof-of-concept implementation

Here's a proof-of-concept implementation, using Python 3 and gmpy2 (configured to calculate Hydra, but it could easily be adapted for other consistent Collatz sequences with $m=2$):

import gmpy2
import time

# The Rules of Hydra:
# if x_n = 2y + 0, then x_{n+1} = 3y + 0, i.e. 2x_{n+1} = 3x_n + 0
# if x_n = 2y + 1, then x_{n+1} = 3y + 1, i.e. 2x_{n+1} = 3x_n - 1
m = 2        # modulus; denominator of the ratio between successive elements
r = 3        # numerator of the ratio between successive elements
J = [0, -1]  # J[x_n % m] is the difference between m*x_{n+1} and r*x_n
x_0 = 3      # first term of the sequence

def f(n):
    """The largest power of 2 that divides into the argument"""
    return n & ~(n - 1)

global mod_m_exp
if m == 2:
    mod_m_exp = lambda v,e: gmpy2.f_mod_2exp(v, e)
else:
    mod_m_exp = lambda v,e: gmpy2.f_mod(v, gmpy2.pow(m, e))

global div_m_exp
if m == 2:
    div_m_exp = lambda v,e: gmpy2.f_div_2exp(v, e)
else:
    div_m_exp = lambda v,e: gmpy2.f_div(v, gmpy2.pow(m, e))

global mul_m_exp
if m == 2:
    mul_m_exp = lambda v,e: gmpy2.mpz(v) << e
else:
    mul_m_exp = lambda v,e: gmpy2.mul(v, gmpy2.pow(m, e))

r_exp_cache = {1: gmpy2.mpz(r)}
def r_exp(e):
    """Returns r raised to the power of e. e must be a power of 2."""
    global r_exp_cache
    if not e in r_exp_cache:
        r_exp_cache[e] = gmpy2.square(r_exp(e/2))
    return r_exp_cache[e]

# The bulk of the calculation is to calculate w_n and j_n for each n.
# The definitions are:
# w_n = x_n mod m**f(n)
# j_n = m**f(n) * x_n - r**f(n) * x_{n-f(n)}
# The output from the program is the sequence of x_n mod m.
# This can be calculated by taking the values of w_n mod m.
w = {0: x_0}  # most recently seen w for each f value
j = {}        # most recently seen j for each f value

modulus_count = {0: 0, 1: 0}

n = 0
perf_counter_timestamp = time.perf_counter_ns()
while True:
    last_f = f(n)
    last_x_mod_m = mod_m_exp(w[last_f], 1)
    # print(last_x_mod_m, end="")

    modulus_count[last_x_mod_m] += 1
    if n % 1000000 == 0:
        last_timestamp = perf_counter_timestamp
        perf_counter_timestamp = time.perf_counter_ns()
        print("Reached n = ", n, "; modulus counts: ", modulus_count,
              "; time for last 1000000 elements = ",
              (perf_counter_timestamp - last_timestamp) // 1000000,
              " ms", sep="")

    n = n + 1
    cur_f = f(n)

    # j can be a sum of multiple terms; J[last_x_mod_m] is always present
    # but if f > 1 there are other terms too
    m_shift = cur_f - 1
    r_shift = 1
    new_j = mul_m_exp(J[last_x_mod_m], m_shift)
    while m_shift >= r_shift:
        m_shift -= r_shift
        new_j += mul_m_exp(j[r_shift] * r_exp(r_shift), m_shift)
        r_shift *= 2
    j[cur_f] = new_j

    # w can be calculated directly from the new j and the appropriate past w
    wrap = lambda v: mod_m_exp(v, cur_f * 2)
    past_w = w[f(n - f(n))]
    new_w = wrap(wrap(new_j) + wrap(wrap(r_exp(cur_f)) * wrap(past_w)))
    assert(mod_m_exp(new_w, cur_f) == 0)
    new_w = div_m_exp(new_w, cur_f)
    w[cur_f] = new_w

This implementation is probably not suitable for serious use due to having poor constant factors: a faster implementation would use an alternative algorithm for low $f(n)$ values. It is also intened only for $m=2$, because the complexity result is dependent on modulus by powers of $m$ being fast regardless of the size of the dividend, but gmpy2 is only able to store numbers in binary and that operation is quick only if the modulus is a power of the base. As such, a full implementation of efficient consistent Collatz calculation would probably involve writing a new arbitrarily-large-integers library which is able to store numbers in arbitrary bases.