A number theory problem where pi appears surprisingly For a given positive integer $M$, the sequence $\{a_n\}$ starts from $a_{2M+1}=M(2M+1)$ and $a_k$ is the largest multiple of $k$ no more than $a_{k+1}+M$, i.e.
$$a_k=k\left\lfloor\frac{a_{k+1}+M}{k}\right\rfloor,\quad k=1,2,\cdots,2M.$$
The original problem asks me to show that $a_1<4M^2$ for $M\ge 3$. Then I write a program to check large $M$s like $M\sim10000000$ and surprisingly find that
$$\lim_{M\rightarrow\infty}\frac{a_1}{M^2}=\pi.$$
Is this true and why does this happen?
 A: Note that $a_k$ is always a multiple of $k$ (in particular, this structure disrupts the random model proposed in comments).  We can exploit this structure to simplify the recurrence and clarify the dynamics by making the change of variables $a_k = k (b_k - M)$ (the shift by $M$ is convenient to eliminate some lower order terms), then $b_{2M+1} = 2M$, the $b_k$ are integers, and for $k=1,\dots,2M$ one has
\begin{align*}
 b_k &= \frac{a_k}{k} + M\\
&= \left\lfloor \frac{a_{k+1}+M}{k} \right\rfloor + M \\
&= \left\lfloor \frac{(k+1)(b_{k+1}-M)+M}{k} \right\rfloor + M\\
&= (b_{k+1}-M) + \left\lfloor \frac{(b_{k+1}-M) + M}{k} \right\rfloor + M\\
&= b_{k+1} + \left\lfloor \frac{b_{k+1}}{k}\right \rfloor
\end{align*}
and so we have arrived at the difference equation
$$ b_{k+1} - b_k = - \left\lfloor \frac{b_{k+1}}{k} \right\rfloor.$$
Thus, as $k$ decrements from $k+1$, $b_k$ will increment by $1$ in the regime $k \leq b_{k+1} < 2k$, by two in the regime $2k \leq b_{k+1} < 3k$, and so forth.  This is quite a stable recurrence and can be analyzed by the standard technique of passing to a rescaled limit and studying the resulting asymptotic differential equation.  For simplicity we argue heuristically.  We expect $b_k$ to have magnitude $M$ when $k \sim M$, so it is natural to introduce a rescaling
$$ b_k = 2M f_M(\frac{k}{2M})$$
then we have the initial condition $f_M(1 + \frac{1}{2M})=1$, and for $k = 2Mt$ for $0 \leq t \leq 1$ a multiple of $\frac{1}{2M}$ we then have after a short calculation
$$ \frac{f_M(t+\frac{1}{2M}) - f_M(t)}{1/2M} = - \left\lfloor \frac{f_M(t+\frac{1}{2M})}{t} \right\rfloor.$$
Formally passing to the limit $M \to \infty$, we thus expect $f_M$ to converge in some suitable sense to a continuous, piecewise differentiable function $f: [0,1] \to {\bf R}$ with $f(1)=1$ that solves the ODE
$$ f'(t) = - \left\lfloor \frac{f(t)}{t} \right\rfloor$$
for all $0 \leq t \leq 1$.
(To make this rigorous, one can for instance rewrite both the difference equation and the presumed limiting ODE as partial sum equations and integral equations respectively to improve the stability, and then apply a compactness theorem such as the Arzela-Ascoli theorem to conclude, in the spirit of constructing various mild or weak solutions in PDE.)  One can intuitively think of the graph of $f$ as the trajectory of a light ray passing through different media that travels in straight lines until it hits one of the flat boundaries $f(t)=t$, $f(t)=2t$, $f(t)=3t$, etc., at which point it "refracts" to a different slope. (A numerical plot of the graphs of $f$ and $f_M$ for some reasonably large $M$ would be quite revealing, but I do not have the time to generate these plots at present.)
One can solve this ODE explicitly from the initial condition $f(1)=1$ as a continuous piecewise linear function.  If we define the sequence $1 = t_0 > t_1 > t_2 > \dots$ by the recursion
$$ t_n = \frac{2n}{2n+1} t_{n-1}$$
so that $t_n$ is a partial half-Wallis product
$$ t_n = \frac{2}{3} \cdot \frac{4}{5} \cdot \dots \cdot \frac{2n}{2n+1} = \frac{(n! 2^n)^2}{(2n+1)!}$$
(which already begins to hint at the emergence of $\pi$) then one can show by induction that
$$ f(t) = n t_{n-1} + n(t_{n-1}-t) = (n+1) t_n - n (t-t_n)$$
for $t_n \leq t \leq t_{n-1}$ and all $n \geq 1$.  In particular
$$ f(t_n) = (n+1) t_n.$$
From Stirling's approximation we have for large $n$ that
$$ t_n = \frac{(n! 2^n)^2}{(2n+1)!} \sim \frac{\sqrt{\pi}}{2\sqrt{n}}.$$
We remark that this asymptotic can also be derived from (and is in fact equivalent to) the Wallis product $\frac{\pi}{2} = (\frac{2}{1} \cdot \frac{2}{3}) \cdot (\frac{4}{3} \cdot \frac{4}{5}) \cdot (\frac{6}{5} \cdot \frac{6}{7}) \cdot \dots$ after some routine algebraic manipulation.
Thus if $k \approx 2M t_n \approx 2M \frac{\sqrt{\pi}}{2\sqrt{n}}$, then
$$ b_k \approx 2M f(t_n) \approx (n+1) 2M \frac{\sqrt{\pi}}{2\sqrt{n}}$$
and hence
\begin{align*}
a_k &\approx k(b_k-M) \\
&\approx 2M \frac{\sqrt{\pi}}{2\sqrt{n}} ((n+1) 2M \frac{\sqrt{\pi}}{2\sqrt{n}} - M )\\
&\approx M^2 (\pi + O(n^{-1/2})).
\end{align*}
Sending $n \to \infty$, we morally obtain
$$ a_1 \approx M^2 \pi$$
as expected.
It is likely that a more careful quantitative accounting of error terms can turn the above heuristic analysis into a rigorous argument, but I leave this to the interested reader.
We remark also that this recurrence can be used to give a (quite inefficient) algorithm for computing $\pi$ on a simple adding machine, basically by implementing the Wallis product formula in a rather pedestrian fashion.
