Is the sequence $(a_n)$ defined by $a_n=\tan{a_{n-1}}$, $a_0=1$, dense in $\Bbb{R}$?
Before I start, I want to point out that I do not answer the specific question. I believe that no one can with the current tools available in mathematics. However, what I will give an answer that is at the state of the art. And here comes the second disclaimer: I will summon quite recent research (about 15 years old), with a wealth of technical definitions and theorems. Giving precise proofs would take a small article, so I will have to be a little sketchy to keep things short. The upside is that I can get quite precise propositions.
For all $x \in \mathbb{R}$, we define recursively $a_n (x)$ by $a_0 (x) = x$ and $a_{n+1} (x) = \tan (x)$. I will not be able to say anything about the sequence $(a_n(1))_{n \geq 0}$, which is why I do not answer your question. Nevertheless, I will be able to say quite a lot of things about $(a_n(x))_{n \geq 0}$ for a Lebesgue-generic real number $x$. Then, you just have to hope that $1$ has these generic properties (there is no obvious reason it has not). This restriction comes from the fact that I work with a chaotic dynamical system, which exhibit sensitivity to initial conditions. For such systems you can usually work out a lot of generic properties, but you cannot tell if any given point is indeed generic or not (and there are indeed a lot of non-generic points and weird potential behaviours).
We can see $(\mathbb{R}, \tan)$ as a dynamical system. It is well-defined up to a countable set, which does not matter in the following since a countable set is negligible for the Lebesgue measure. The first thing to note is that the $\tan$ function is $\pi$-periodic. Hence, to dynamical system goes to the quotient. Let $I := (-\pi/2, \pi/2)$, let $\pi : \ \mathbb{R} \to I$ be the reduction modulo $\pi$, and let:
\begin{equation*} T : \ \left\{ \begin{array}{lll} I & \to & I \\ x & \mapsto & \pi \circ \tan(x) \end{array} \right. . \end{equation*}
Then, for all $n \geq 0$,
\begin{equation*} a_{n+1} = \tan \circ T^n \circ \pi. \end{equation*}
Hence, the density of $(a_n(x))$ in $\mathbb{R}$ is equivalent to the density of $T^n (x)$ in $I$. We can go further, and relate the equidistribution of $T^n (x)$ in $I$ with the equidistribution of $(a_n(x))$ in $\mathbb{R}$ (more on that later).
Now, we only need to study the system $(I, T)$. What does it look like? Here is its plot (courtesy of Wolfram Alpha):
This system has a few nice properties, and a few not-so-nice properties. It has countably many branches (not so nice), but these branches are surjective, which ensure a Markov property and a big image property (nice). It is piecewise $\mathcal{C}^2$, which ensures a bounded distortion property. Hence, the fact that it has countably many branches is manageable. At first sight, it share a lot of properties with the Gauss map $x \to \{1/x\}$, which has been thoroughly studied:
Now, $T$ has a real not-so-nice property. We like our maps to be expanding, that is, that $\inf |T'| > 1$: that gets us positive Lyapunov exponents, positive entropy, and chaos in any sense of the word. The Gauss map is expanding everywhere except at $1$, where its derivative is $-1$. But that is not a hard problem; you can check that if you iterate the Gauss map, what you get is expanding. For $T$, however, the problem is more serious: the derivative at $0$ is $1$, but $0$ is a fixed point. Hence, $(T^n)' (0) = 1$ for all $n$. A point which lands close to $0$ will be eventually repelled (that's because $|\tan (x) > x|$ for non-zero $x$), but that may take a lot of time. Actually, the orbit spend most of their time close to zero. This is an expanding map with a neutral fixed point, and it belongs to a class of dynamical systems which has been studied at length in the past twenty years.
The phenomena which arise are about the same as for one of the Liverani-Saussol-Vaienti maps [LSV] with a parameter $\alpha = 2$:
\begin{equation*} T_2 : \ \left\{ \begin{array}{lll} (0,1] & \to & (0,1] \\ x & \mapsto & \left\{ \begin{array}{lll} x(1+4x^2) & \text{ if } & x \in (0,1/2],\\ 2x-1 & \text{ if } & x \in (1/2,1] \end{array}\right. \end{array} \right. . \end{equation*}
The choice of the parameter $\alpha = 2$ is due to the fact that $T_2 (x) = x + 4x^3 + o(x^3)$ at $0$, just as $T (x) = x + x^3/3 + o(x^3)$ at $0$: the second-order terms are of the same order.
One way of studying these LSV maps, or $T$ is to induce the system away from the neutral fixed point (here $0$). For instance, let $A := (-\pi/2, -\pi/4) \cup (\pi/4, \pi/2)$, where $-\pi/4$ and $\pi/4$ are chosen as the end-points of the branch of $T$ containing $0$. Let $\varphi(x) := \inf_{n > 0} \{ T^n (x) \in A\}$, and let $T_A (x) := T^{\varphi (x)} (x)$. By looking closely at what happens around $0$ (points are repelled), it can be shown that $\varphi (x) < + \infty$ for all but countably many $x$. Hence, $T_A$ is well-defined almost everywhere. Finally, $T_A$ is expanding. Taking into account the Markov property, the big image property and the bounded distortion property, $(A, T_A)$ can be endowed with a Gibbs-Markov structure (see [A], Chapter 4, and in particular Sections 4.7 for definitions and 4.8 for an application to interval maps).
We can then show that $T_A$ has a unique absolutely continuous invariant probability measure, with respect to which it is ergodic and mixing (literature says has finitely many ergodic components of positive Lebesgue measure, and the surjectivity of the branches of $T$ gives transitivity and aperiodicity, hence mixing), and which has full support. In particular, for almost every $x \in A$, the orbit $(T_A^n (x))_{n \geq 0}$ is equidistributed (for this invariant measure), and thus dense.
The map $(I, T)$ can then be seen as a Rokhlin tower over $(A, T_A)$ of height $\varphi$. This is already sufficient to prove that almost every $x \in I$ has dense orbits in $I$, and thus that $(a_n (x))_{n \geq 0}$ is dense in $\mathbb{R}$ for Lebesgue-almost every $x$. But we can go further!
Using [LSV], we can show that the transformation $(I, T)$ has, up to positive multiplicative constant, a unique positive $\sigma$-finite measure $\mu$ which is $T$-invariant and absolutely continuous with respect to the Lebesgue measure. Its density has a piecewise continuous version $g$, which is bounded away from $0$, and such that $g(x) = \Theta (x^{-2})$ at $0$. The dynamical system $(I, T, \mu)$ is ergodic. The measure $\mu$ is isomorphic the the natural invariant measure on the Rokhlin tower over $A$; but the height of the tower has infinite mean (because the trajectories spend a lot of time to get away from $0$), so the measure $\mu$ is infinite.
Then, $\nu := \tan_* \mu$ is a $\sigma$-finite measure on $\mathbb{R}$ which is $\tan$-invariant. It also has a piecewise continuous density $h$, such that $h(x) = \Theta (x^{-2})$ at $0$ and $h(x) \sim c x^{-2}$ at $\pm \infty$ for some $c > 0$, and which satisfies the functional equation:
\begin{equation*} \sum_{n \in \mathbb{Z}} h (\arctan (x) + n \pi) = (1+x^2) h(x). \end{equation*}
To finish with this answer, here are a few consequences.
Density: for almost every $x \in \mathbb{R}$, the sequence $(a_n (x))_{n \geq 0}$ is dense in $\mathbb{R}$.
Trajectories spend most of their time close to $0$: for almost every $x \in \mathbb{R}$, for any neighborhood $B$ of $0$,
\begin{equation*} \lim_{n \to + \infty} \frac{\# \{0 \leq k < n : \ a_k (x) \in B\}}{n} = 1. \end{equation*}
- Local time: there is a constant $C$ such that for all measurable $B$ bounded away from $0$, for all probability measure $\mathbb{P}$ absolutely continuous with respect to the Lebesgue measure on $\mathbb{R}$,
\begin{equation*} \lim_{n \to + \infty} \frac{\# \{0 \leq k < n : \ a_k \in B\}}{\sqrt{n}} = C \int_B h(t) \ dt \cdot |\mathcal{N}|, \end{equation*}
where $\mathcal{N}$ is a Gaussian random variable of variance $\pi/2$ and the convergence is in distribution in $(\mathbb{R}, \mathbb{P})$. This comes from [A], Corollary 3.7.3. This tells you roughly that about $1/\sqrt{n}$ of the $(a_k (x))_{0 \leq k < n}$ will be in any subset bounded away from zero, with a significant stochastic variation. That can make for nice numerical experiments, I you wish so: compute $\# \{0 \leq k < n : \ |a_k| (x) > N \}/\sqrt{n}$ for fixed large $n$ and many random values of $x$ (so as to get a statistical approximation of its distribution), and look at what happens when you let $N$ grow.
- It is likely (to my knowledge, this specific result is not in the literature, but I do not think it is that hard to show) that a typical trajectory $(\# \{0 \leq k < N : \ a_k (x) \in B\})_{n \geq 0}$ will behave at large scale like a typical trajectory of the local time of a $1$-dimensional simple random walk.
We can hope that such properties are true for $(a_n (1))_{n \geq 0}$.
[A] J. Aaronson, An introduction to infinite ergodic theory, Mathematical Surveys and Monographs, Vol.50, American Mathematical Society, 1997.
[LSV] C. Liverani, B. Saussol and S. Vaienti, A probabilistic approach to intermittency, Ergodic Theory and Dynamical Systems 19 (1999), 671--685, available here.