Hybrid Numerical Integration Algorithms

A hybrid single precision algorithm can use the following strategy: use (3) for $\vert\rho\vert \le 0.8$ and use (6) for $\vert\rho\vert \ge 0.8$. Extensive testing shows that the maximum absolute error for this hybrid algorithm is $2.5\cdot10^{-7}$, if a 5-point Gauss-Legendre integration rule (see Davis and Rabinowitz, 1984, pp. 95-101) is used. A quadruple precision BVN algorithm based on Owen's (1956) series method was used as a source of highly accuracy BVN values for all of the tests described in this section.

The present author investigated the use of more Gauss points to produce an algorithm that would have maximum errors near $10^{-16}$ (double precision). In order to avoid very high numbers of Gauss points when $\vert\rho\vert \simeq 1$, a higher order Taylor expansion can be used in (5). Adding one more term to this expansion yields

$$\frac{1}{\sqrt{1-x^2}} e^{-\frac{shk}{1+\sqrt{1-x^2}}}= e^{-\frac{shk}{2}}(1+(4-shk)x^2(1 + (12-shk)x^2/16)/8) +O(x^6).$$ (7)

If this modification is incorporated into equation (6) then a 20-point Gauss-Legendre rule used with the modified equation (6) for $\vert\rho\vert \geq 0.925$, and used with equation (3) for $\vert\rho\vert < 0.925$ results in a maximum absolute error less than $5\cdot10^{-16}$ in double precision. A more efficient algorithm can be constructed by selecting Gauss rules with fewer points for the smaller values of $\vert\rho\vert$. It was found, for example, that if the algorithm is modified to use only six points for $\vert\rho\vert < 0.3$ and twelve points for $0.3 \leq \vert\rho\vert < 0.75$ then the same level of accuracy could be maintained. These final modifications were incorporated into the Fortran implementation that is available from the author. Further refinements are obviously possible, at the expense of additional implementation complexity.