nag_elliptic_integral_rf (s21bbc) calculates an approximation to the integral
where
,
,
and at most one is zero.
The basic algorithm, which is due to Carlson (1978) and Carlson (1988), is to reduce the arguments recursively towards their mean by the rule:
remaining third intermediate value argument. (This ordering, which is possible because of the symmetry of the function, is done for technical reasons related to the avoidance of overflow and underflow.)
and the function may be approximated adequately by a 5th-order power series:
where
.
The truncation error involved in using this approximation is bounded by
and the recursive process is stopped when this truncation error is negligible compared with the machine precision.
Within the domain of definition, the function value is itself representable for all representable values of its arguments. However, for values of the arguments near the extremes the above algorithm must be modified so as to avoid causing underflows or overflows in intermediate steps. In extreme regions arguments are pre-scaled away from the extremes and compensating scaling of the result is done before returning to the calling program.