...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/special_functions/ellint_rf.hpp>
namespace boost { namespace math { template <class T1, class T2, class T3> calculatedresulttype ellint_rf(T1 x, T2 y, T3 z) template <class T1, class T2, class T3, class Policy> calculatedresulttype ellint_rf(T1 x, T2 y, T3 z, const Policy&) }} // namespaces
#include <boost/math/special_functions/ellint_rd.hpp>
namespace boost { namespace math { template <class T1, class T2, class T3> calculatedresulttype ellint_rd(T1 x, T2 y, T3 z) template <class T1, class T2, class T3, class Policy> calculatedresulttype ellint_rd(T1 x, T2 y, T3 z, const Policy&) }} // namespaces
#include <boost/math/special_functions/ellint_rj.hpp>
namespace boost { namespace math { template <class T1, class T2, class T3, class T4> calculatedresulttype ellint_rj(T1 x, T2 y, T3 z, T4 p) template <class T1, class T2, class T3, class T4, class Policy> calculatedresulttype ellint_rj(T1 x, T2 y, T3 z, T4 p, const Policy&) }} // namespaces
#include <boost/math/special_functions/ellint_rc.hpp>
namespace boost { namespace math { template <class T1, class T2> calculatedresulttype ellint_rc(T1 x, T2 y) template <class T1, class T2, class Policy> calculatedresulttype ellint_rc(T1 x, T2 y, const Policy&) }} // namespaces
#include <boost/math/special_functions/ellint_rg.hpp>
namespace boost { namespace math { template <class T1, class T2, class T3> calculatedresulttype ellint_rg(T1 x, T2 y, T3 z) template <class T1, class T2, class T3, class Policy> calculatedresulttype ellint_rg(T1 x, T2 y, T3 z, const Policy&) }} // namespaces
These functions return Carlson's symmetrical elliptic integrals, the functions have complicated behavior over all their possible domains, but the following graph gives an idea of their behavior:
The return type of these functions is computed using the result type calculation rules when the arguments are of different types: otherwise the return is the same type as the arguments.
template <class T1, class T2, class T3> calculatedresulttype ellint_rf(T1 x, T2 y, T3 z) template <class T1, class T2, class T3, class Policy> calculatedresulttype ellint_rf(T1 x, T2 y, T3 z, const Policy&)
Returns Carlson's Elliptic Integral R_{F}:
Requires that all of the arguments are nonnegative, and at most one may be zero. Otherwise returns the result of domain_error.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
template <class T1, class T2, class T3> calculatedresulttype ellint_rd(T1 x, T2 y, T3 z) template <class T1, class T2, class T3, class Policy> calculatedresulttype ellint_rd(T1 x, T2 y, T3 z, const Policy&)
Returns Carlson's elliptic integral R_{D}:
Requires that x and y are nonnegative, with at most one of them zero, and that z >= 0. Otherwise returns the result of domain_error.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
template <class T1, class T2, class T3, class T4> calculatedresulttype ellint_rj(T1 x, T2 y, T3 z, T4 p) template <class T1, class T2, class T3, class T4, class Policy> calculatedresulttype ellint_rj(T1 x, T2 y, T3 z, T4 p, const Policy&)
Returns Carlson's elliptic integral R_{J}:
Requires that x, y and z are nonnegative, with at most one of them zero, and that p != 0. Otherwise returns the result of domain_error.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
When p < 0 the function returns the Cauchy principal value using the relation:
template <class T1, class T2> calculatedresulttype ellint_rc(T1 x, T2 y) template <class T1, class T2, class Policy> calculatedresulttype ellint_rc(T1 x, T2 y, const Policy&)
Returns Carlson's elliptic integral R_{C}:
Requires that x > 0 and that y != 0. Otherwise returns the result of domain_error.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
When y < 0 the function returns the Cauchy principal value using the relation:
template <class T1, class T2, class T3> calculatedresulttype ellint_rg(T1 x, T2 y, T3 z) template <class T1, class T2, class T3, class Policy> calculatedresulttype ellint_rg(T1 x, T2 y, T3 z, const Policy&)
Returns Carlson's elliptic integral R_{G}:
Requires that x and y are nonnegative, otherwise returns the result of domain_error.
The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.
There are two sets of tests.
Spot tests compare selected values with test data given in:
B. C. Carlson, Numerical computation of real or complex elliptic integrals. Numerical Algorithms, Volume 10, Number 1 / March, 1995, pp 1326.
Random test data generated using NTL::RR at 1000bit precision and our implementation checks for roundingerrors and/or regressions.
There are also sanity checks that use the interrelations between the integrals to verify their correctness: see the above Carlson paper for details.
These functions are computed using only basic arithmetic operations, so there isn't much variation in accuracy over differing platforms. Note that only results for the widest floatingpoint type on the system are given as narrower types have effectively zero error. All values are relative errors in units of epsilon.
Table 8.58. Error rates for ellint_rc
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Microsoft Visual C++ version 14.1 


RC: Random data 
Max = 0ε (Mean = 0ε) 
Max = 0.995ε (Mean = 0.433ε) 
Max = 0.962ε (Mean = 0.407ε) 
Table 8.59. Error rates for ellint_rd
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Microsoft Visual C++ version 14.1 


RD: Random data 
Max = 0ε (Mean = 0ε) 
Max = 2.73ε (Mean = 0.831ε) 
Max = 2.16ε (Mean = 0.803ε) 
RD: y = z 
Max = 0.896ε (Mean = 0.022ε) 
Max = 2.65ε (Mean = 0.82ε) 
Max = 16.5ε (Mean = 0.843ε) 
RD: x = y 
Max = 0.824ε (Mean = 0.0272ε) 
Max = 2.85ε (Mean = 0.865ε) 
Max = 3.51ε (Mean = 0.816ε) 
RD: x = 0, y = z 
Max = 0ε (Mean = 0ε) 
Max = 1.19ε (Mean = 0.522ε) 
Max = 1.16ε (Mean = 0.497ε) 
RD: x = y = z 
Max = 0ε (Mean = 0ε) 
Max = 0.998ε (Mean = 0.387ε) 
Max = 1.03ε (Mean = 0.418ε) 
RD: x = 0 
Max = 0ε (Mean = 0ε) 
Max = 2.79ε (Mean = 0.883ε) 
Max = 2.64ε (Mean = 0.894ε) 
Table 8.60. Error rates for ellint_rg
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Microsoft Visual C++ version 14.1 


RG: Random Data 
Max = 0.983ε (Mean = 0.0172ε) 
Max = 3.95ε (Mean = 0.951ε) 
Max = 3.65ε (Mean = 0.929ε) 
RG: two values 0 
Max = 0ε (Mean = 0ε) 
Max = 0ε (Mean = 0ε) 
Max = 0ε (Mean = 0ε) 
RG: All values the same or zero 
Max = 0ε (Mean = 0ε) 
Max = 0.992ε (Mean = 0.288ε) 
Max = 1.06ε (Mean = 0.348ε) 
RG: two values the same 
Max = 0.594ε (Mean = 0.0103ε) 
Max = 1.51ε (Mean = 0.404ε) 
Max = 1.96ε (Mean = 0.374ε) 
RG: one value zero 
Max = 0ε (Mean = 0ε) 
Max = 2.14ε (Mean = 0.722ε) 
Max = 1.96ε (Mean = 0.674ε) 
Table 8.61. Error rates for ellint_rf
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Microsoft Visual C++ version 14.1 


RF: Random data 
Max = 0ε (Mean = 0ε) 
Max = 2.54ε (Mean = 0.674ε) 
Max = 2.02ε (Mean = 0.677ε) 
RF: x = y = z 
Max = 0ε (Mean = 0ε) 
Max = 0.991ε (Mean = 0.345ε) 
Max = 0.999ε (Mean = 0.34ε) 
RF: x = y or y = z or x = z 
Max = 0.536ε (Mean = 0.00658ε) 
Max = 1.95ε (Mean = 0.418ε) 
Max = 1.21ε (Mean = 0.394ε) 
RF: x = 0, y = z 
Max = 0ε (Mean = 0ε) 
Max = 0.894ε (Mean = 0.338ε) 
Max = 0.999ε (Mean = 0.407ε) 
RF: z = 0 
Max = 0ε (Mean = 0ε) 
Max = 1.7ε (Mean = 0.539ε) 
Max = 1.89ε (Mean = 0.587ε) 
Table 8.62. Error rates for ellint_rj
GNU C++ version 7.1.0 
GNU C++ version 7.1.0 
Microsoft Visual C++ version 14.1 


RJ: Random data 
Max = 0.52ε (Mean = 0.0184ε) 
Max = 186ε (Mean = 6.67ε) 
Max = 215ε (Mean = 7.66ε) 
RJ: 4 Equal Values 
Max = 0ε (Mean = 0ε) 
Max = 0.998ε (Mean = 0.387ε) 
Max = 1.03ε (Mean = 0.418ε) 
RJ: 3 Equal Values 
Max = 0ε (Mean = 0ε) 
Max = 20.8ε (Mean = 0.986ε) 
Max = 39.9ε (Mean = 1.17ε) 
RJ: 2 Equal Values 
Max = 0.6ε (Mean = 0.0228ε) 
Max = 220ε (Mean = 6.64ε) 
Max = 214ε (Mean = 5.28ε) 
RJ: Equal z and p 
Max = 0.742ε (Mean = 0.0166ε) 
Max = 17.2ε (Mean = 1.16ε) 
Max = 16.1ε (Mean = 1.14ε) 
The key of Carlson's algorithm [Carlson79] is the duplication theorem:
By applying it repeatedly, x, y, z get closer and closer. When they are nearly equal, the special case equation
is used. More specifically, [R F] is evaluated from a Taylor series expansion to the fifth order. The calculations of the other three integrals are analogous, except for R_{C} which can be computed from elementary functions.
For p < 0 in R_{J}(x, y, z, p) and y < 0 in R_{C}(x, y), the integrals are singular and their Cauchy principal values are returned via the relations: