Ray tracing plays a key role in lens design area, and is an important tool to research problems in physics like optics. Nowadays, ray tracing becomes ubiquitous and is widely used in optical automatic design discipline, such as aberration analysis, optimization, and tolerance calculation. Due to impulse of application requirements, optical systems like space camera develop towards directions for large scale, high degree of accuracy and complication. The magnitude of aberrations increases exponentially as growths on focal length and apertures, even a minor perturbation error can result in severe degeneration of image quality. As a consequence, the stringent requirements of precision, accuracy and stability of ray tracing are getting higher. Reliable commercial software, for example, America’s Zemax, has high precision in ray tracing, because of commercial purpose, the process of ray tracing is a black box. It is now more important to understand what error factors are formed for ray tracing, and how to effectively reduce these running errors. In this paper, from floating point arithmetic perspective, error model for ray tracing is provided. This error model is not only suited for meridional rays, but also skew rays. Starting from IEEE Standard for Binary Floating-Point Arithmetic, presentation error and rounding error are analyzed, followed by the computation process of ray’s intersection point with a quadratic surface, then rounding error expression for the intersection point is presented. In addition, error expression for distance along the ray from the reference surface to the next surface is also induced. These two error expressions are called error model, and it clearly indicates that spatial coordinates on the reference surface, direction vector and distance between the two adjacent surfaces are the main error sources. Based on the error model, some of effective measures, for instance, reprojection, spatial transformation, and direction vector’s normalization are taken to reduce the rounding error. Moreover, during the process of solving quadratic equation, in order to avoid substantial increase in relative error called catastrophic cancellation, conjugate number method is utilized. Numerical experiments and classical optical design for space camera are also given. From numerical computing view, two MPFR (multiple-precision floating-point reliable library) based precision tests are introduced to verify our method mathematically, the experimental results show that our algorithm has the same precision (14 significant digits) as MPFR, while the existing method fails to pass tests, and only has 8 significant digits at most. Moreover, both the Cassegrain space camera and off-axis three-mirror-anastigmat space camera are used to illustrate our method’s accuracy, experimental results indicate that our method has higher precision, which has more than 5 to 6 orders of magnitudes than existing method. In addition, our algorithm has higher precision than commercial optical design software Zemax, and residuals are less than 3 orders of magnitudes on average than Zemax.