Ray tracing plays a key role in lens design area, and it is an important tool to study the problems in physics like optics. Nowadays, ray tracing becomes ubiquitous and is widely used in optical automatic design, such as aberration analysis, optimization, and tolerance calculation. With the impulse of application requirements, optical systems like space camera develop towards large scale, high degree of accuracy and complication. The magnitude of aberrations increases exponentially with the growth of focal length and aperture, even a minor perturbation error can result in severe degeneration of image quality. As a consequence, the stringent requirements for precision, accuracy and stability of ray tracing turn 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 these running errors can be reduced effectively. In this paper, from floating point arithmetic perspective, an error model for ray tracing is provided. This error model is suitable for not only 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, in the process of solving quadratic equation, conjugate number method is utilized in order to avoid increasing substantially in relative error called catastrophic cancellation. Numerical experiments and classical optical design for space camera are also given. From numerical computing view, two precision tests based on Multiple Precision Floating-Point Reliable (MPFR) library 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 has only 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, more than 5 to 6 orders of magnitudes than the existing method. In addition, our algorithm has higher precision than the commercial optical design software Zemax, and residuals are 3 orders of magnitudes on average less than Zemax.