  In this paper, we are concerned with symmetric integrators for the nonlinear
relativistic Klein--Gordon (NRKG) equation with a dimensionless parameter
$0<\varepsilon\ll 1$, which is inversely proportional to the speed of light.
The highly oscillatory property in time of this model corresponds to the
parameter $\varepsilon$ and the equation has strong nonlinearity when $\eps$ is
small. There two aspects bring significantly numerical burdens in designing
numerical methods. We propose and analyze a novel class of symmetric
integrators which is based on some formulation approaches to the problem,
Fourier pseudo-spectral method and exponential integrators. Two practical
integrators up to order four are constructed by using the proposed symmetric
property and stiff order conditions of implicit exponential integrators. The
convergence of the obtained integrators is rigorously studied, and it is shown
that the accuracy in time is improved to be $\mathcal{O}(\varepsilon^{3}
\hh^2)$ and $\mathcal{O}(\varepsilon^{4} \hh^4)$ for the time stepsize $\hh$.
The near energy conservation over long times is established for the multi-stage
integrators by using modulated Fourier expansions. These theoretical results
are achievable even if large stepsizes are utilized in the schemes. Numerical
results on a NRKG equation show that the proposed integrators have improved
uniform error bounds, excellent long time energy conservation and competitive
efficiency.
