Boost C++ Libraries

...one of the most highly regarded and expertly designed C++ library projects in the world. Herb Sutter and Andrei Alexandrescu, C++ Coding Standards

PrevUpHomeNext

Elliptic Integrals - Carlson Form

Synopsis
#include <boost/math/special_functions/ellint_rf.hpp>
namespace boost { namespace math {

template <class T1, class T2, class T3>
calculated-result-type ellint_rf(T1 x, T2 y, T3 z)

template <class T1, class T2, class T3, class Policy>
calculated-result-type 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>
calculated-result-type ellint_rd(T1 x, T2 y, T3 z)

template <class T1, class T2, class T3, class Policy>
calculated-result-type 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>
calculated-result-type ellint_rj(T1 x, T2 y, T3 z, T4 p)

template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type 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>
calculated-result-type ellint_rc(T1 x, T2 y)

template <class T1, class T2, class Policy>
calculated-result-type 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>
calculated-result-type ellint_rg(T1 x, T2 y, T3 z)

template <class T1, class T2, class T3, class Policy>
calculated-result-type ellint_rg(T1 x, T2 y, T3 z, const Policy&)

}} // namespaces
Description

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>
calculated-result-type ellint_rf(T1 x, T2 y, T3 z)

template <class T1, class T2, class T3, class Policy>
calculated-result-type ellint_rf(T1 x, T2 y, T3 z, const Policy&)

Returns Carlson's Elliptic Integral RF:

Requires that all of the arguments are non-negative, 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>
calculated-result-type ellint_rd(T1 x, T2 y, T3 z)

template <class T1, class T2, class T3, class Policy>
calculated-result-type ellint_rd(T1 x, T2 y, T3 z, const Policy&)

Returns Carlson's elliptic integral RD:

Requires that x and y are non-negative, 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>
calculated-result-type ellint_rj(T1 x, T2 y, T3 z, T4 p)

template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ellint_rj(T1 x, T2 y, T3 z, T4 p, const Policy&)

Returns Carlson's elliptic integral RJ:

Requires that x, y and z are non-negative, 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>
calculated-result-type ellint_rc(T1 x, T2 y)

template <class T1, class T2, class Policy>
calculated-result-type ellint_rc(T1 x, T2 y, const Policy&)

Returns Carlson's elliptic integral RC:

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>
calculated-result-type ellint_rg(T1 x, T2 y, T3 z)

template <class T1, class T2, class T3, class Policy>
calculated-result-type ellint_rg(T1 x, T2 y, T3 z, const Policy&)

Returns Carlson's elliptic integral RG:

Requires that x and y are non-negative, 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.

Testing

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 13-26.

Random test data generated using NTL::RR at 1000-bit precision and our implementation checks for rounding-errors and/or regressions.

There are also sanity checks that use the inter-relations between the integrals to verify their correctness: see the above Carlson paper for details.

Accuracy

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 floating-point type on the system are given as narrower types have effectively zero error. All values are relative errors in units of epsilon.

Table 6.58. Error rates for ellint_rc

Microsoft Visual C++ version 12.0
Win32
double

GNU C++ version 5.1.0
linux
double

GNU C++ version 5.1.0
linux
long double

Sun compiler version 0x5130
Sun Solaris
long double

RC: Random data

Max = 0.962ε (Mean = 0.407ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 2.4ε (Mean = 0.624ε))

Max = 0.995ε (Mean = 0.433ε)

Max = 0.995ε (Mean = 0.438ε)


Table 6.59. Error rates for ellint_rd

Microsoft Visual C++ version 12.0
Win32
double

GNU C++ version 5.1.0
linux
double

GNU C++ version 5.1.0
linux
long double

Sun compiler version 0x5130
Sun Solaris
long double

RD: Random data

Max = 2.16ε (Mean = 0.803ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 2.59ε (Mean = 0.878ε))

Max = 2.73ε (Mean = 0.831ε)

Max = 2.73ε (Mean = 0.829ε)

RD: y = z

Max = 16.5ε (Mean = 0.843ε)

Max = 0.896ε (Mean = 0.022ε)

(GSL 1.16: Max = 2.88ε (Mean = 0.839ε))

Max = 2.65ε (Mean = 0.82ε)

Max = 2.65ε (Mean = 0.819ε)

RD: x = y

Max = 3.51ε (Mean = 0.816ε)

Max = 0.824ε (Mean = 0.0272ε)

(GSL 1.16: Max = 3.74ε (Mean = 0.84ε))

Max = 2.85ε (Mean = 0.865ε)

Max = 2.85ε (Mean = 0.865ε)

RD: x = 0, y = z

Max = 1.16ε (Mean = 0.493ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 2ε (Mean = 0.656ε))

Max = 1.19ε (Mean = 0.522ε)

Max = 1.19ε (Mean = 0.522ε)

RD: x = y = z

Max = 1.03ε (Mean = 0.418ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 1.03ε (Mean = 0.418ε))

Max = 0.998ε (Mean = 0.387ε)

Max = 0.998ε (Mean = 0.387ε)

RD: x = 0

Max = 2.64ε (Mean = 0.894ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 2.85ε (Mean = 0.781ε))

Max = 2.79ε (Mean = 0.883ε)

Max = 2.79ε (Mean = 0.883ε)


Table 6.60. Error rates for ellint_rg

Microsoft Visual C++ version 12.0
Win32
double

GNU C++ version 5.1.0
linux
double

GNU C++ version 5.1.0
linux
long double

Sun compiler version 0x5130
Sun Solaris
long double

RG: Random Data

Max = 3.65ε (Mean = 0.929ε)

Max = 0.983ε (Mean = 0.0172ε)

(GSL 1.16: Max = 0.983ε (Mean = 0.0172ε))

Max = 3.95ε (Mean = 0.951ε)

Max = 3.95ε (Mean = 0.951ε)

RG: two values 0

Max = 0ε (Mean = 0ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 0ε (Mean = 0ε))

Max = 0ε (Mean = 0ε)

Max = 0ε (Mean = 0ε)

RG: All values the same or zero

Max = 1.06ε (Mean = 0.348ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 0ε (Mean = 0ε))

Max = 0.992ε (Mean = 0.288ε)

Max = 0.992ε (Mean = 0.288ε)

RG: two values the same

Max = 1.96ε (Mean = 0.374ε)

Max = 0.594ε (Mean = 0.0103ε)

(GSL 1.16: Max = 0.594ε (Mean = 0.0103ε))

Max = 1.51ε (Mean = 0.404ε)

Max = 1.51ε (Mean = 0.404ε)

RG: one value zero

Max = 1.96ε (Mean = 0.674ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 0ε (Mean = 0ε))

Max = 2.14ε (Mean = 0.722ε)

Max = 2.14ε (Mean = 0.722ε)


Table 6.61. Error rates for ellint_rf

Microsoft Visual C++ version 12.0
Win32
double

GNU C++ version 5.1.0
linux
double

GNU C++ version 5.1.0
linux
long double

Sun compiler version 0x5130
Sun Solaris
long double

RF: Random data

Max = 2.02ε (Mean = 0.677ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 2.73ε (Mean = 0.804ε))

Max = 2.54ε (Mean = 0.674ε)

Max = 2.54ε (Mean = 0.674ε)

RF: x = y = z

Max = 0.999ε (Mean = 0.335ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 0.999ε (Mean = 0.34ε))

Max = 0.991ε (Mean = 0.345ε)

Max = 0.991ε (Mean = 0.345ε)

RF: x = y or y = z or x = z

Max = 1.21ε (Mean = 0.394ε)

Max = 0.536ε (Mean = 0.00658ε)

(GSL 1.16: Max = 2.89ε (Mean = 0.749ε))

Max = 1.95ε (Mean = 0.418ε)

Max = 1.57ε (Mean = 0.418ε)

RF: x = 0, y = z

Max = 0.999ε (Mean = 0.407ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 1.29ε (Mean = 0.527ε))

Max = 0.894ε (Mean = 0.338ε)

Max = 0.894ε (Mean = 0.338ε)

RF: z = 0

Max = 1.89ε (Mean = 0.587ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 2.54ε (Mean = 0.781ε))

Max = 1.7ε (Mean = 0.539ε)

Max = 1.7ε (Mean = 0.539ε)


Table 6.62. Error rates for ellint_rj

Microsoft Visual C++ version 12.0
Win32
double

GNU C++ version 5.1.0
linux
double

GNU C++ version 5.1.0
linux
long double

Sun compiler version 0x5130
Sun Solaris
long double

RJ: Random data

Max = 119ε (Mean = 4.32ε)

Max = 0.52ε (Mean = 0.0184ε)

(GSL 1.16: Max = 3.57ε (Mean = 0.704ε) And other failures.)

Max = 186ε (Mean = 6.67ε)

Max = 186ε (Mean = 6.7ε)

RJ: 4 Equal Values

Max = 1.03ε (Mean = 0.418ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 1.03ε (Mean = 0.418ε))

Max = 0.998ε (Mean = 0.387ε)

Max = 0.998ε (Mean = 0.387ε)

RJ: 3 Equal Values

Max = 39.9ε (Mean = 1.12ε)

Max = 0ε (Mean = 0ε)

(GSL 1.16: Max = 3.96ε (Mean = 1.06ε))

Max = 20.8ε (Mean = 0.986ε)

Max = 18.2ε (Mean = 0.917ε)

RJ: 2 Equal Values

Max = 214ε (Mean = 5.05ε)

Max = 0.6ε (Mean = 0.0228ε)

(GSL 1.16: Max = 2.57ε (Mean = 0.754ε))

Max = 220ε (Mean = 6.64ε)

Max = 135ε (Mean = 5.3ε)

RJ: Equal z and p

Max = 15.4ε (Mean = 1.05ε)

Max = 0.742ε (Mean = 0.0166ε)

(GSL 1.16: Max = 2.62ε (Mean = 0.699ε))

Max = 17.2ε (Mean = 1.16ε)

Max = 16.6ε (Mean = 1.15ε)


Implementation

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 RC which can be computed from elementary functions.

For p < 0 in RJ(x, y, z, p) and y < 0 in RC(x, y), the integrals are singular and their Cauchy principal values are returned via the relations:


PrevUpHomeNext