6
$\begingroup$

Sequence A76132 starting as $1,1,2,4,10,36,218,\ldots$ of the OEIS is recursively defined by $a(1)=1$ and $a(n)=\sum_{k=1}^{n-1}a(n-k)^k$ for $n\geq 2$.

It is eventually periodic of period 1,1 and 34 modulo $2,3$ and $5$ (the proof is essentially trivial: it is enough to check eventual periodicity on a big enough finite initial chunk).

Curiously, this seems to stop: either the non-periodic part is very long (of length at least 100) or the period is of length longer than 5000 for $p=7$.

Moreover, the sequence seems to be eventually periodic modulo $2^k$: up to $2^3=8$ the sequence is eventually constant modulo $2^k$ (for $k=1,2,3$) after that it seems to be eventually periodic with period-length $2^{k-3}$.

Similarly, it seems to be eventually periodic with period length $3^k$ for $k\geq 2$.

Curiously, it does not seem to be eventually periodic modulo $25$.

Is there an easy explanation for this behaviour?

Added: The generating series $\sum a(n)t^n$ is algebraic modulo every prime number. Eventual periodicity modulo $p$ boils down to showing that it is a rational series modulo $p$.

Final addition. The remark of reuns implies that $\sum a(n)t^n$ is always rational modulo any integer $m$. The aperiodic part and the period length are however generally of order $O(m^m)$ which explains why they are not easily observable.

$\endgroup$
2
  • 5
    $\begingroup$ Pick some integer $m$, let $b(n)_j=\sum_{k=1}^n (a_{n-k})^{k+j} \bmod m$ for $j\in 0\ldots 2\phi(m)$. If for some $n_1\ne n_2$, $b(n_1)=b(n_2)$ then $(a_n)_{n\ge n_1}$ is $n_2-n_1$-periodic modulo $m$. ​Of course such $n_1,n_2$ must exist as $b(n)$ takes its values in a finite set $\endgroup$
    – reuns
    Jan 6, 2022 at 20:47
  • $\begingroup$ Great observation! I think one can diminish the dimension of the vectors $b(n)$ somewhat. This implies however that the aperiodic and the periodic part are huge in general. $\endgroup$ Jan 6, 2022 at 21:15

1 Answer 1

2
$\begingroup$

The comment from "reuns" is the real answer-- the sequence is periodic for all $M$. For prime $M$, we can slightly improve the implied bound as follows.

We track the pair: $$\left(j, \sum_{k=1}^{n}(a_{n-k}){k+j}\right) \bmod M$$ for $j\in 0,1,..., \phi(M)-1$ (instead of $2\phi(M)$). We may end up detecting a multiple of a cycle, but the total size of the state space (and thus maximum cycle length) is smaller: $\phi(M)M^{\phi(M)}$.

This is still not a great bound, e.g., the bound for M=11 is 259,374,246,010. In practice, though, the periods may be shorter. I'm including a little snippet of Python code for finding these. Running time is a fraction of a second.

First, the results:

Extending series to 1000000 places
###   M=2   ###
  minimal period: start 1, length 1
###   M=3   ###
  minimal period: start 5, length 1
###   M=5   ###
  minimal period: start 3, length 37
###   M=7   ###
  minimal period: start 253, length 807
###   M=11   ###
  minimal period: start 303832, length 9
###   M=13   ###
  ... Period not found ...
###   M=17   ###
  ... Period not found ...
###   M=19   ###
  ... Period not found ...

Next, the code:

# Try all the values from 2 to 20
modulus_list = range(2,21)

# "max_cycle_length" is optional.  If set to 0, we detect cycles of
# any length, but memory grows linearly as the experiment runs.  If
# it's >0, then memory is bounded and we can extend the cycle
# indefinitely given enough time.
max_cycle_length = 0

import numpy as np
from math import gcd

def phi(n):
    amount = 0        
    for k in range(1, n + 1):
        if gcd(n, k) == 1:
            amount += 1
    return amount

# Encode the historical state into a single, hashable number
def encode(ar):
    out=0
    multiplier=1
    for a in ar:
        out+=a*multiplier
        multiplier*=M
    return((out, step % phiM))

num_trials = 1000000
print('Extending series to %d places' % num_trials)
if max_cycle_length>0:
    print('Maximum detectable cycle length: %d' % max_cycle_length)
for M in modulus_list:
    phiM = phi(M)
    # Restrict to primes for now
    if phiM < M-1:
        continue
    print("###   M=%d   ###" % M)

    state=np.ones(phiM)
    history={}
    a=np.zeros(num_trials)

    periodic = False
    for step in range(num_trials):
        new_value = state[step % phiM]
        a[step] = new_value
        powered = new_value
        for i in range(phiM):
            j = (1+i+step) % phiM
            state[j] = (powered+state[j]) % M
            powered = (powered*new_value) % M
        enc = encode(state)
        if enc in history:
            periodic = True
            break
        else:
            history[enc] = step
        if max_cycle_length > 0:
            if len(history) == max_cycle_length:
                history = {}
    if not periodic:
        print('  ... Period not found ...')
        continue

    # We may have a multiple of the fundamental period; minimize it.
    start = history[enc]+1
    end = step
    if end <= start+1:
        print('  minimal period: start %d, length 1' % (start))
    else:
        chunk = a[start:start+28]
        for i in range(start+2, end-len(chunk)+1):
            if np.array_equal(chunk, a[i:i+len(chunk)]):
                print('  minimal period: start %d, length %d' % (start, i))
                break
        else:
            print('  minimal period: start %d, length %d' % (start, end-start))

Extending num_trials=1000000 did not turn up anything else, but a little more patience might reveal something

$\endgroup$
1
  • $\begingroup$ This is the sense of my 'Final addition'. Thanks for the additional data: $\endgroup$ Jan 7, 2022 at 20:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.