I don't know any VASP-specific details, but the basic answer is that the three versions are almost identical, and solve the same equations in the same way. If we take vasp_std as the reference, then the differences are:
vasp_gam
If your system has time-reversal symmetry (true for most simulations), then $+k$ and $-k$ are symmetry-related and we can restrict our calculation to just one half of the Brillouin Zone. At the $\Gamma$-point, $k=(0,0,0)$ and so $+k=-k$; in this case, time-reversal symmetry relates $+G$ to $-G$, where $G$ is a reciprocal lattice vector. This means that the plane-wave coefficients $c(G)$ and $c(-G)$ are related by symmetry (in fact, they are the complex conjugate of each other). Thus, if we're using the $\Gamma$-point we only need to consider half of the plane-waves explicitly, and we can generate the rest by symmetry. This halves the memory requirements of the calculation, and makes the calculation go twice as quickly.
One consequence of this symmetry between $+G$ and $-G$ is that the wavefunctions in real-space can be made to have zero imaginary components -- i.e. entirely real. This means that any inner product or Hermitian matrix element must also be real, since we could calculate them in real-space if we wanted and now everything is real there. Thus, whenever we calculate a matrix element we only need to bother computing the real part, since we already know that the imaginary part will be zero; this reduces the time to compute matrix elements by an additional factor of two.
Perhaps the best part, is that the cubically scaling parts of the calculation, such as orthogonalisation of the Kohn-Sham wavefunction, now involves cubically scaling steps on a real matrix, whereas before they were complex. These operations are now $2^3=8$ times faster, meaning that the asymptotic speed-up when using the $\Gamma$-point version is x8.
NB A slightly more subtle, but very similar symmetry argument can be made for points at the boundary of the Brillouin Zone ($X$-points etc) and some DFT software can also exploit that, e.g. ABINIT.
vasp_ncl
Spin-orbit coupling is a relativistic effect which is included in the external potential in the (electronic) Kohn-Sham equations, and hence enters via the pseudopotential term in VASP and similar programs (e.g. ABINIT, CASTEP or Quantum Espresso). The pseudopotential generator has to solve an atomic Dirac-like equation now, not a Schroedinger-like one, and this means that the solutions are spinors. A full solution needs 4-component spinors, but in fact for ground state calculations the major components are so large compared to the minor ones, that 2-component spinors suffice for the valence states. The Kohn-Sham wavefunctions thus have twice the plane-wave coefficients as in a "normal" calculation, and require twice the memory storage.
The usual approach to handling the spin-orbit interaction is from dal Corso and Conte, and moves all the spin-orbit terms into the non-local pseudopotential projectors; this has the consequence that the local pseudopotential term must not be for a chemically-relevant angular-momentum channel. The pseudopotential projectors now depend not on the orbital angular momentum $l$ directly, but on the total angular momentum $j=l+s$. This requires there to be at least twice as many pseudopotential projectors for a spin-orbit pseudopotential, as for a "normal" pseudopotential.
The Kohn-Sham valence wavefunctions are 2-component spinors, so they gain a spinor index in addition to the plane-wave, band and k-point indices (NB there is no $s_z$ index, since $s_z$ is no longer a good quantum number and spin is entirely described by the spinors). Inner products now have to be calculated by summing over spinor indices as well, e.g.,
\begin{align}\tag{1}
\langle \psi_{bk}\vert\phi_{b^\prime k}\rangle =& \sum_{\sigma} \int \psi_{b\sigma k}^\ast(r)\phi_{b^\prime\sigma k}(r)d^3 r.
\end{align}
Note that we only needed one spinor index here, $\sigma$, because the spinor basis is orthogonal and we only wanted the inner product.
For matrix elements the situation is more complicated, because some interactions do couple different spinor components, and so they gain two spinor indices, $\sigma$ and $\sigma^\prime$. Perhaps the simplest example is actually not a matrix element exactly, but the spin-density-matrix $\rho$, which becomes,
\begin{align}\tag{2}
\rho_{\sigma \sigma^\prime}(r)=\sum_{bk} \psi_{b\sigma k}^\ast(r)\psi_{b\sigma^\prime k}(r).
\end{align}
The trace of this matrix over the spinor components is the usual charge density, but the other elements describe the spin, which may now point in any direction in real-space -- i.e. it is a vector. (You can calculate the components of spin in any particular direction using the $2\times 2$ Pauli spin matrices acting on this density matrix.)
The addition of two spinor indices to several key objects in the Kohn-Sham scheme means that for large calculations the computational cost is 8 times that of a "normal" calculation (4 times compared to a collinear-spin calculation).
(Note that in the presence of spin-orbit coupling, the Kohn-Sham states are complex in real-space even at the $\Gamma$-point, unless the spin is aligned along $z$ everywhere.)
