Sequence A76132 eventually periodic modulo $2,3$ and $5$

Question

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.

Answer 1

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