I will show how to generally evaluate $a_1^{\,a_2^{\,\cdots}} \bmod m$, recursively.
If $(a_1 \bmod m) \leq 1$, we have $a_1^{\,a_2^{\,\cdots}} \equiv a_1 \mod m$.
Otherwise, if $\gcd(a_1, m) = 1$, we have $a_1^{\,a_2^{\,\cdots}} \equiv a_1^{\,a_2^{\,\cdots} \bmod \phi(m)}\mod m$. Recursively evaluate $r = a_2^{\,a_3^{\,\cdots}} \bmod \phi(m)$, and then calculate $a_1^r \bmod m$.
Otherwise, factorize $m$ into primes. If $m$ is a prime power $p^k$, but $\gcd(a_1, p^k) \neq 1$, then $\gcd(a_1, p^k) = p^i$. Write $a'_1 = a_1/p^i$, then we have $a_1^{\,a_2^{\,\cdots}} \equiv a_1'^{\,a_2^{\,\cdots}}\cdot(p^i)^{\,a_2^{\,\cdots}} \mod p^k$.
We recursively evaluate $r \equiv a_1'^{\,a_2^{\,\cdots}} \equiv a_1'^{\,a_2^{\,\cdots}\bmod \phi(p^k)} \equiv a_1'^{\,a_2^{\,\cdots}\bmod p^k - p^{k-1}} \mod p^k$.
$s \equiv (p^i)^{\,a_2^{\,\cdots}} \equiv p^{i\,a_2^{\,\cdots}} \bmod p^k$ is different from the rest. If $i\,a_2^{\,a_3^{\,\cdots}} \geq k$ then $s = 0$, otherwise we simply need to evaluate bigger and bigger prefixes of $a$. If either the power tower gets too big, or we reach a $0$ or $1$, we stop, evaluate the power tower and find $s$. This is always a finite process, the tower will not be longer than $1 + \log_2k$ (this is even a bad upper bound, the real one is the tetra-logarithm).
Then we find $st \bmod m$ and we're done.
Recursively evaluate $a_1^{\,a_2^{\,\cdots}} \bmod p^k$ for each $p,k$ in the factorization of $m$. Then use the Chinese remainder theorem to find $a_1^{\,a_2^{\,\cdots}} \bmod m$.
Note that for any sequence $a$, and any modulus $m$, this terminates, as at every recursive step $m$ decreases, and $m$ is finite. Even when $a$ is infinite.
An implementation in Python using sympy (using zero based indexing for $a$):
from math import gcd, log
from sympy.ntheory import totient, factorint
from sympy.ntheory.modular import crt
def evaltowermax(l, k):
r = 1
for e in l[::-1]:
# Prevent evaluation of large powers, 0.1 to account for errors.
if log(e)*r - 0.1 > log(k):
r = k
break
r = e**r
return r
def modulartower(af, m, n=0):
a = af(n); g = gcd(a, m)
if a % m <= 1: return a % m
if g == 1: return pow(a, modulartower(af, totient(m), n + 1), m)
f = factorint(m)
if len(f) == 1: # Prime power.
p, i = factorint(g).popitem()
k = f[p]
tower = [af(ti) for ti in range(n, n + k.bit_length() + 1)]
s = pow(p, i * evaltowermax(tower, k), m)
if s == 0: return 0
aprimef = lambda l: a // p**i if l == n else af(l)
t = modulartower(aprimef, p**k - p**(k-1), n)
return s*t % m
m = [p**k for p, k in f.items()]
r = [modulartower(af, p**k, n) for p, k in f.items()]
return crt(m, r)[0]
print(modulartower(lambda n: 2017 + n, 10**20))
This computes the last 20 digits of $2017^{2018^{\cdots}}$ in an instant as $77345043177395978241$.
Simplified algorithm due to user Feersum:
If $(a_1 \bmod m) \leq 1$, we have $a_1^{\,a_2^{\,\cdots}} \equiv a_1 \mod m$.
Otherwise, factor $m$ into primes $p_i$. We calculate $\displaystyle x = \prod_{\gcd(p_i, a) = 1} p_i$ and $y = m/x$.
Then, we compute $a_1^{\,a_2^{\,\cdots}}\bmod x$ and $a_1^{\,a_2^{\,\cdots}}\bmod y$ and use the Chinese remainder theorem to find $a_1^{\,a_2^{\,\cdots}}\bmod xy = a_1^{\,a_2^{\,\cdots}}\bmod m$.
Since $a_1$ and $x$ are coprime we have $a_1^{\,a_2^{\,\cdots}}\equiv a_1^{\,a_2^{\,\cdots} \bmod \phi(x)} \mod x$. We recursively compute $r = a_2^{\,a_3^{\,\cdots}} \bmod \phi(x)$ and then compute $a_1^r \mod x$ directly.
$a_1$ and $y$ is a bit interesting, as $y$ only consists of primes that are found in the decomposition of $a_1$. So for a large enough $k$ we have $a_1^k \equiv 0 \mod y$. For an infinite sequence $a$ without $0$ or $1$ elements, this is always the case. If an infinite sequence contains $0$ or $1$ elements or the sequence is finite, we must evaluate the prime tower $a_1^{\,a_2^{\,\cdots}}\bmod y$, luckily we only need to evaluate $\log_2 y + 1$ steps at worst.
