What is the difference between controllability and reachability?

They are not the same problem, but equivalent for linear continuous systems. Also, reachability and controllability gramians are slightly different. To understand the differences let's start with a general linear continuous time-varying system. $$ \dot{x}=A(t)x(t)+B(t)u(t) $$

Its solution can be given as $$ x(t) = \phi(t,t_0) x(t_0) + \int_{t_0}^t \phi(t,\tau) B(\tau) u(\tau) d\tau $$ where $\phi(\cdot,\cdot)$ is the state transition matrix.

Now, let's say we want to "reach" to the state $x(t_1)=x_1$ at time $t_1$ for a given $x(t_0)=x_0$. Then, we can use the input function $$ u(t) = B^T(t) \phi^T(t_1, t) W_r^{-1}(t_1,t_0) \left(x_1 - \phi(t_1,t_0) x_0 \right) $$ where $$W_r(t_1,t_0) := \int_{t_0}^{t_1} \phi(t_1, \eta) B(\eta) B^T(\eta) \phi^T(t_1, \eta) d\eta$$

Note that if the reachability gramian is full-rank, we can reach any state we want from any initial condition, hence full reachability. If it does not have a full-rank you can still show that the reachable subspace at time $t_1$ is $$\begin{align*}\mathcal{R}(t_0;t_1) &= \operatorname{Im} \int_{t_0}^{t_1} \phi(t_1,\tau) B(\tau) d\tau \\ &= \operatorname{Im} W_r(t_1,t_0)\end{align*}$$

For controllability, suppose your final state is given as $x(t_1)=x_1$ and you want to find which initial states can reach this final state. Then using the properties of the state transition matrix, $$x_0 = \phi^{-1}(t_1,t_0) x_1 - \int_{t_0}^{t_1} \phi(t_0, \tau) B(\tau) u(\tau) d\tau$$ which is now essentially the same problem with reachability, but backwards in time. So, the controllable subspace is $$\begin{align*}\mathcal{C}(t_0;t_1) &= \operatorname{Im} \int_{t_0}^{t_1} \phi(t_0,\tau) B(\tau) d\tau \\ &= \operatorname{Im} W_c(t_1,t_0)\end{align*}$$ where $$W_c(t_1,t_0) := \int_{t_0}^{t_1} \phi(t_0, \eta) B(\eta) B^T(\eta) \phi^T(t_0, \eta) d\eta$$

Discrete time case is more interesting, because reachability and controllability is not equivalent in this case as you pointed out. The reason is the state transition matrix (which is $A^k$ for discrete LTI case) might not be invertible (we cannot always go backwards in time) as it is in continuous time case. But the thought process is the same.

To summarize,

• For full reachability in linear continuous systems: $\operatorname{Im} W_r(t_1,t_0) = \mathbb{R}^n$
• For full controllability in linear continuous systems: $\mathbb{R}^n = \operatorname{Im} \phi(t_1,t_0) \subseteq \operatorname{Im} W_c(t_1,t_0)$
• For full reachability in LTI discrete systems: $\operatorname{Im} \sum_{i=0}^{k-1} A^i B = \mathbb{R}^n$
• For full controllability in LTI discrete systems: $\operatorname{Im}A^n \subseteq \operatorname{Im} \sum_{i=0}^{n-1} A^i B$