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.

We present a symbolic-numeric Las Vegas algorithm for factoring Fuchsian ordinary differential operators with rational function coefficients. The new algorithm combines ideas of van Hoeij's "local-to-global" method and of the "analytic" approach proposed by van der Hoeven. It essentially reduces to the former in "easy" cases where the local-to-global method succeeds, and to an optimized variant of the latter in the "hardest" cases, while handling intermediate cases more efficiently than both.

Let S\subset \bR^n be a compact basic semi-algebraic set defined as the real solution set of multivariate polynomial inequalities with rational coefficients. We design an algorithm which takes as input a polynomial system defining S and an integer p\geq 0 and returns the n-dimensional volume of S at absolute precision 2^-p . Our algorithm relies on the relationship between volumes of semi-algebraic sets and periods of rational integrals. It makes use of algorithms computing the Picard-Fuchs differential equation of appropriate periods, properties of critical points, and high-precision numerical integration of differential equations. The algorithm runs in essentially linear time with respect to~p. This improves upon the previous exponential bounds obtained by Monte-Carlo or moment-based methods. Assuming a conjecture of Dimca, the arithmetic cost of the algebraic subroutines for computing Picard-Fuchs equations and critical points is singly exponential in n and polynomial in the maximum degree of the input.

We introduce an algorithm to compare a binary floating-point (FP) number and a decimal FP number, assuming the “binary encoding” of the decimal formats is used, and with a special emphasis on the basic interchange formats specified by the IEEE 754-2008 standard for FP arithmetic. It is a two-step algorithm: a first pass, based on the exponents only, quickly eliminates most cases, then, when the first pass does not suffice, a more accurate second pass is performed. We provide an implementation of several variants of our algorithm, and compare them.

This work introduces an approach to the computer-assisted implementation of mathematical functions geared toward special functions such as those occurring in mathematical physics. The general idea is to start with an exact symbolic representation of a function and automate as much as possible of the process of implementing it. In order to deal with a large class of special functions, our symbolic representation is an implicit one: the input is a linear differential equation with polynomial coefficients along with initial values. The output is a C program to evaluate the solution of the equation using domain splitting, argument reduction and polynomial approximations in double-precision arithmetic, in the usual style of mathematical libraries. Our generation method combines symbolic-numeric manipulations of linear ODEs with interval-based tools for the floating-point implementation of "black-box" functions. We describe a prototype code generator that can automatically produce implementations on moderately large intervals. Implementations on the whole real line are possible in some cases but require manual tool setup and code integration. Due to this limitation and as some heuristics remain, we refer to our method as "semi-automatic" at this stage. Along with other examples, we present an implementation of the Voigt profile with fixed parameters that may be of independent interest.

We introduce a software-oriented algorithm that allows one to quickly compare a binary64 floating-point (FP) number and a decimal64 FP number, assuming the "binary encoding" of the decimal formats specified by the IEEE 754-2008 standard for FP arithmetic is used. It is a two-step algorithm: a first pass, based on the exponents only, makes it possible to quickly eliminate most cases, then when the first pass does not suffice, a more accurate second pass is required. We provide an implementation of several variants of our algorithm, and compare them.

The series expansion at the origin of the Airy function Ai(x) is alternating and hence problematic to evaluate for x > 0 due to cancellation. Based on a method recently proposed by Gawronski, Müller, and Rein hard, we exhibit two functions F and G, both with nonnegative Taylor expansions at the origin, such that Ai(x) = G(x)/F(x). The sums are now well-conditioned, but the Taylor coefficients of G turn out to obey an ill-conditioned three-term recurrence. We use the classical Miller algorithm to overcome this issue. We bound all errors and our implementation allows an arbitrary and certified accuracy, that can be used, e.g., for providing correct rounding in arbitrary precision.

We present a new algorithm for computing hyperexponential solutions of linear ordinary differential equations with polynomial coefficients. The algorithm relies on interpreting formal series solutions at the singular points as analytic functions and evaluating them numerically at some common ordinary point. The numerical data is used to determine a small number of combinations of the formal series that may give rise to hyperexponential solutions.

This article describes the implementation in the software package NumGfun of classical algorithms that operate on solutions of linear differential equations or recurrence relations with polynomial coefficients, including what seems to be the first general implementation of the fast high-precision numerical evaluation algorithms of Chudnovsky & Chudnovsky. In some cases, our descriptions contain improvements over existing algorithms. We also provide references to relevant ideas not currently used in NumGfun.