Quadratic forms in Gaussian random variables yield generalized noncentral chi-square distributions which appear in many test statistics. Classical cumulative distribution function (cdf) evaluation methods (e.g., characteristic-function inversion, saddlepoint approximations) can be effective but are not designed for reliable finite-precision computation and offer limited complexity insight. We study finite-precision cdf evaluation and its bit complexity.

The finite precision implementation of mathematical functions frequently depends on polynomial approximations. A key characteristic of this approach is that rounding errors occur both when representing the coefficients of the polynomial on a finite number of bits, and when evaluating it in finite precision arithmetic. Hence, to find a best polynomial, for a given fixed degree, norm, and interval, it is necessary to account for both the approximation error and the floating-point evaluation error. While efficient algorithms were already developed for taking into account the approximation error, the evaluation part is usually a posteriori handled, in an ad hoc manner. Here, we formulate a semi-infinite linear optimization problem whose solution is a best polynomial with respect to the supremum norm of the sum of both errors. This problem is then solved with an iterative exchange algorithm, which can be seen as an extension of the well-known Remez exchange algorithm. An open source C implementation using the Sollya library is presented and tested on several examples, which are then analyzed and compared against state-of-the-art Sollya routines.

We present an error analysis of an algorithm due to Serra et al. (Journal of Guidance Control and Dynamics, 2016) for computing the orbital collision probability in the short term encounter model. The algorithm reduces the numerical computation of the collision probability to that of the sum of a series whose coefficients are produced by a linear recurrence relation, and is specifically designed to avoid cancellation issues in the evaluation of the sum. While its numerical stability was observed experimentally and a truncation error bound was derived, the evaluation error was not studied. Here we give a rigorous bound on the accumulated rounding error when Serra et al.’s algorithm is implemented in floating-point arithmetic. For a unit roundoff u and a truncation order N, the bound is of the form (N + A)u + o(u), where o(u) stands for small compared to u, explicitly bounded terms. The constant A explicitly depends on the problem parameters and dominates N in practice. Our analysis is based on the observation that the generating series of the errors affecting each individual term is solution to a perturbed form of a differential equation satisfied by the Laplace transform of a function related to the collision probability.

Machine implementation of mathematical functions often relies on polynomial approximations. The particularity is that rounding errors occur both when representing the polynomial coefficients on a finite number of bits, and when evaluating it in finite precision. Hence, for finding the best polynomial (for a given fixed degree, norm and interval), one has to consider both types of errors: approximation and evaluation. While efficient algorithms were already developed for taking into account the approximation error, the evaluation part is usually a posteriori handled, in an ad-hoc manner. Here, we formulate a semi-infinite linear optimization problem whose solution is the best polynomial with respect to the supremum norm of the sum of both errors. This problem is then solved with an iterative exchange algorithm, which can be seen as an extension of the well-known Remez algorithm. A discussion and comparison of the obtained results on different examples are finally presented.