We present several numerical methods and establish their error estimates for
the discretization of the nonlinear Dirac equation in the nonrelativistic limit
regime, involving a small dimensionless parameter $0<\varepsilon\ll 1$ which is
inversely proportional to the speed of light. In this limit regime, the
solution is highly oscillatory in time, i.e. there are propagating waves with
wavelength $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. We
begin with the conservative Crank-Nicolson finite difference (CNFD) method and
establish rigorously its error estimate which depends explicitly on the mesh
size $h$ and time step $\tau$ as well as the small parameter $0<\varepsilon\le
1$. Based on the error bound, in order to obtain `correct' numerical solutions
in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the CNFD method
requests the $\varepsilon$-scalability: $\tau=O(\varepsilon^3)$ and
$h=O(\sqrt{\varepsilon})$. Then we propose and analyze two numerical methods
for the discretization of the nonlinear Dirac equation by using the Fourier
spectral discretization for spatial derivatives combined with the exponential
wave integrator and time-splitting technique for temporal derivatives,
respectively. Rigorous error bounds for the two numerical methods show that
their $\varepsilon$-scalability is improved to $\tau=O(\varepsilon^2)$ and
$h=O(1)$ when $0<\varepsilon\ll 1$ compared with the CNFD method. Extensive
numerical results are reported to confirm our error estimates.