To your first question: No, it is false. For $n = 256 - 13$, the number $a_n$ has $\nu_2\left(a_n\right) = 7 < 8 = \nu_2 \left(n+13\right)$.
HOWEVER, it is almost correct: namely, it is correct whenever $128 \nmid n+13$ (so the smallest counterexample is $n = 128 - 13 = 115$). This is the following theorem:
Theorem 1. We have $\nu_2\left(a_n\right) = \nu_2\left(n+13\right)$ for all $n \in \mathbb{N}$ that don't satisfy $128 \mid n+13$.
Proof of Theorem 1 (sketched). We observe the following facts:
Each nonnegative integer $n$ and each positive integer $k$ satisfy $a_{n+k} \equiv a_n \mod k$. (This is proven by induction on $n$. The base case boils down to $a_k \equiv 1 \mod k$, which follows from $a_k = k a_{k-1} + 1$. The induction step uses the recursion $a_n = na_{n-1} + 1$.)
Each $n \in \mathbb{N}$ satisfies $a_{n+128} \equiv a_n \mod 128$. (This follows by applying fact 1 to $k = 128$.)
Each $n \in \mathbb{N}$ satisfies $128 \nmid a_n$ unless $128 \mid n+13$. (This needs only to be checked for the first $128$ values of $n$, due to fact 2 above.)
Combining facts 2 and 3, obtain $\nu_2\left(a_{n+128}\right) = \nu_2\left(a_n\right)$ for each $n \in \mathbb{N}$ unless $128 \mid n+13$.
Using fact 4, the Theorem 1 can be proven by induction (with $128$ base cases, unfortunately). $\blacksquare$
So what about the cases when $128 \mid n+13$ ? Their behavior is weird. The smallest such case is $n = 115$, in which $\nu_2 \left(a_n\right) = 9 > 7 = \nu_2 \left( n + 13 \right)$. However, for all larger $n$'s of the form $2^k - 13$, the sign flips:
Theorem 2. We have $\nu_2\left(a_n\right) = 7 < \nu_2\left(n+13\right)$ for all $n$ of the form $n = 2^k - 13$ with $k \geq 8$.
Proof of Theorem 2 (sketched). We claim that $a_{2^k - 13} \equiv 2^7 \mod 2^8$ for each $k \geq 8$. This is proven by induction on $k$. The induction base ($k = 8$) is verified by computer; the induction step relies on fact 1 from the proof of Theorem 1. $\blacksquare$
Theorem 2 should not suggest that $\nu_2\left(a_n\right)$ is always at most $7$; it is not. Higher values of $\nu_2\left(a_n\right)$ appear when $n+13$ is a multiple of $128$ but not itself a power of $2$. For example, $\nu_2\left(a_{2675}\right) = 15$.
Note that the same inductive argument that we used to prove fact 1 can be used to show the following stronger result:
Proposition 3. For any nonnegative integer $n$ and any positive integer $k$, we have
\begin{align*}
a_{n+k}-a_{n} & \equiv%
\begin{cases}
k, & \text{if }n+k\text{ or }n\text{ is odd};\\
0, & \text{if neither }n+k\text{ nor }n\text{ is odd}%
\end{cases}
\mod 2k.
\end{align*}
This might come useful.
Code: It is not immediately obvious how to compute $a_n$ for large values of $n$ efficiently. After all, $a_n \geq n!$, so we'd at least need a large integer type. Fortunately, all we care about are the $2$-adic valuations $\nu_2\left(a_n\right)$ and the remainders of $a_n$ modulo powers of $2$. The former valuations can easily be obtained from the latter remainders, so we only need to care about the remainders. And we can easily compute the remainders of $a_n$ modulo any given integer $k$ by recursion, because the recursive relation $a_n = na_{n-1} + 1$ descends to $\mathbb{Z}/k\mathbb{Z}$. Here, for example, is some simple Sage (or Python2) code that recursively $a_n$ modulo $2^k$ for given $n$ and $k$:
def a(n, k): # This gives `a_n` modulo `2^k`.
if n == 0:
return 1
return (n * a(n-1) + 1) % (2 ** k)
This is already quite useful (e.g., you can let it compute "a(256-13, 8)" to check that $\nu_2\left(a_n\right) = 7$ for $n = 256-13$), but not really optimal, because it's recursive and eventually (around $n = 1000$) exceeds Python's maximum recursion depth. So let us replace the recursion by dynamic programming.
Here is some better Sage (or Python2) code, which tabulates $\nu_2(a_n)$ for $n = 1, 2, \ldots, m$ given a positive integer $m$:
def val2(k):
# Return the `2`-adic valuation `\nu_2(k)` of the integer `k`.
if k == 0:
return # returns ``None``, standing for `-\infty`.
res = 0
kk = k
while kk % 2 == 0:
res += 1
kk = kk // 2
return res
def aas(n, k):
# Return the list `[a_0, a_1, \ldots, a_{n-1}]` modulo `2^k`.
pwk = 2 ** k
res = [1] * n
for i in range(1, n):
res[i] = (i * res[i-1] + 1) % pwk
return res
def output(m):
# Return a LaTeX table of `\nu_2(a_n)` for `n = 1, 2, \ldots, m`.
k = 3
# Starting with modulo `2^3`, will later replace by better `2`-adic precision.
nums = aas(m+1, k)
while (0 in nums):
k += 1
nums = aas(m+1, k)
vs = [val2(i) for i in nums]
res = r"\begin{array}{|c|c|} \hline" + "\n" + r"n & \nu_2(a_n) \\ \hline "
for i in range(1, m+1):
res += "\n "
res += str(i) + r"&" + str(vs[i]) + r" \\"
res += " \n" + r"\hline \end{array}"
return res
If I let it execute "print output(25)", it returns me the LaTeX code for a table of values of $\nu_2(a_n)$... which I'm not quoting here, because Theorem 1 renders it obsolete. A better idea is to tabulate only those values that Theorem 1 does not cover:
def output2(m):
# Return a LaTeX table of `\nu_2(a_{128n - 13})` for `n = 1, 2, \ldots, m`.
k = 3
# Starting with modulo `2^3`, will later replace by better `2`-adic precision.
M = 128*m - 13
nums = aas(M+1, k)
while (0 in nums):
k += 1
nums = aas(M+1, k)
vs = [val2(i) for i in nums]
res = r"\begin{array}{|c|c|c|} \hline" + "\n" + r"n & \nu_2(a_{128n-13}) & \nu_2(128n) \\ \hline "
for j in range(1, m+1):
i = 128 * j - 13
res += "\n "
res += str(j) + r"&" + str(vs[i]) + r"&" + str(val2(128 * j)) + r" \\"
res += " \n" + r"\hline \end{array}"
return res
Now, executing "print output2(25)" yields
\begin{equation}
\begin{array}{|c|c|c|} \hline
n & \nu_2(a_{128n-13}) & \nu_2(128n) \\ \hline
1&9&7 \\
2&7&8 \\
3&8&7 \\
4&7&9 \\
5&11&7 \\
6&7&8 \\
7&8&7 \\
8&7&10 \\
9&9&7 \\
10&7&8 \\
11&8&7 \\
12&7&9 \\
13&10&7 \\
14&7&8 \\
15&8&7 \\
16&7&11 \\
17&9&7 \\
18&7&8 \\
19&8&7 \\
20&7&9 \\
21&15&7 \\
22&7&8 \\
23&8&7 \\
24&7&10 \\
25&9&7 \\
\hline \end{array} .
\end{equation}
