Fast Switch and Spline Function Inversion Algorithm with Multistep Optimization and k-Vector Search for Solving Kepler’s Equation in Celestial Mechanics

Obtaining the inverse of a nonlinear monotonic function f(x) over a given interval is a common problem in pure and applied mathematics, the most famous example being Kepler’s description of orbital motion in the two-body approximation. In traditional numerical approaches, this problem is reduced to solving the nonlinear equation f(x)−y=0 in each point y of the co-domain. However, modern applications of orbital mechanics for Kepler’s equation, especially in many-body problems, require highly optimized numerical performance. Ongoing efforts continually attempt to improve such performance. Recently, we introduced a novel method for computing the inverse of a one-dimensional function, called the fast switch and spline inversion (FSSI) algorithm. It works by obtaining an accurate interpolation of the inverse function f−1(y) over an entire interval with a very small generation time. Here, we describe two significant improvements with respect to the performance of the original algorithm. First, the indices of the intervals for building the spline are obtained by k-vector search combined with bisection, thereby making the generation time even smaller. Second, in the case of Kepler’s equation, a multistep method for the optimized calculation of the breakpoints of the spline polynomial was designed and implemented in Cython. We demonstrate results that accurately solve Kepler’s equation for any value of the eccentricity e∈[0,1−ϵ], with ϵ=2.22×10−16, which is the limiting error in double precision. Even with modest current hardware, the CPU generation time for obtaining the solution with high accuracy in a large number of points of the co-domain can be kept to around a few nanoseconds per point.

Solving Kepler’s equation with CORDIC double iterations

ABSTRACT In previous work, we developed the idea to solve Kepler’s equation with a CORDIC-like algorithm, which does not require any division, but still requires multiplications in each iteration. Here we overcome this major shortcoming and solve Kepler’s equation using only bitshifts, additions and one initial multiplication. We prescale the initial vector with the eccentricity and the scale correction factor. The rotation direction is decided without correction for the changing scale. We find that double CORDIC iterations are self-correcting and compensate for possible wrong rotations in subsequent iterations. The algorithm needs 75 per cent more iterations and delivers the eccentric anomaly and its sine and cosine terms times the eccentricity. The algorithm can also be adopted for the hyperbolic case. The new shift-and-add algorithm brings Kepler’s equation close to hardware and allows it to be solved with cheap and simple hardware components.

CORDIC-like method for solving Kepler’s equation

Context. Many algorithms to solve Kepler’s equations require the evaluation of trigonometric or root functions. Aims. We present an algorithm to compute the eccentric anomaly and even its cosine and sine terms without usage of other transcendental functions at run-time. With slight modifications it is also applicable for the hyperbolic case. Methods. Based on the idea of CORDIC, our method requires only additions and multiplications and a short table. The table is independent of eccentricity and can be hardcoded. Its length depends on the desired precision. Results. The code is short. The convergence is linear for all mean anomalies and eccentricities e (including e = 1). As a stand-alone algorithm, single and double precision is obtained with 29 and 55 iterations, respectively. Half or two-thirds of the iterations can be saved in combination with Newton’s or Halley’s method at the cost of one division.