A numerical algorithm for solving full gyro orbit of relativistic charged particle motion in magnetized plasmas is presented. The algorithm developed here achieves the following features simultaneously. (1) The time advancement is explicit, and (2) the integration is performed with respect to the observation time in a laboratory frame. (3) It is suitable for accurate long time integration with its volume preserving properties in the phase space. (4) The algorithm can properly treat the E x B drift velocity in electromagnetic fields for large relativistic factors, and (5) can be extended to arbitrary high orders with the aid of symmetric composition methods. Because our algorithm is formulated in the Lorentz covariant form, explicit conservation of the Minkowski norm is not assumed. Nevertheless, no secular growth of the numerical errors in the norm does occur, and its influence can be minimized up to the levels of round-off errors when a high-order scheme is applied. Numerical results are compared with explicit and implicit Runge-Kutta methods. The numerical accuracy and the computational efficiency are discussed for long time integration in toroidal magnetic field configuration. (C) 2017 Elsevier B.V. All rights reserved.