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 Integral D - Legendre Form

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

template <class T1, class T2>
calculated-result-type ellint_d(T1 k, T2 phi);

template <class T1, class T2, class Policy>
calculated-result-type ellint_d(T1 k, T2 phi, const Policy&);

template <class T1>
calculated-result-type ellint_d(T1 k);

template <class T1, class Policy>
calculated-result-type ellint_d(T1 k, const Policy&);

}} // namespaces
Description

These two functions evaluate the incomplete elliptic integral D(φ, k) and its complete counterpart D(k) = D(π/2, k).

The return type of these functions is computed using the result type calculation rules when the arguments are of different types: when they are the same type then the result is the same type as the arguments.

template <class T1, class T2>
calculated-result-type ellint_d(T1 k, T2 phi);

template <class T1, class T2, class Policy>
calculated-result-type ellint_3(T1 k, T2 phi, const Policy&);

Returns the incomplete elliptic integral:

Requires k2sin2(phi) < 1, otherwise returns the result of domain_error (outside this range the result would be complex).

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>
calculated-result-type ellint_d(T1 k);

template <class T1, class Policy>
calculated-result-type ellint_d(T1 k, const Policy&);

Returns the complete elliptic integral D(k) = D(π/2, k)

Requires -1 <= k <= 1 otherwise returns the result of domain_error (outside this range the result would be complex).

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.

Accuracy

These functions are trivially computed in terms of other elliptic integrals and generally have very low error rates (a few epsilon) unless parameter φ is very large, in which case the usual trigonometric function argument-reduction issues apply.

Table 8.66. Error rates for ellint_d (complete)

GNU C++ version 7.1.0
linux
double

GNU C++ version 7.1.0
linux
long double

Sun compiler version 0x5150
Sun Solaris
long double

Microsoft Visual C++ version 14.1
Win32
double

Elliptic Integral E: Mathworld Data

Max = 0.637ε (Mean = 0.368ε)

Max = 1.27ε (Mean = 0.735ε)

Max = 1.27ε (Mean = 0.735ε)

Max = 0.637ε (Mean = 0.368ε)

Elliptic Integral D: Random Data

Max = 0ε (Mean = 0ε)

Max = 1.27ε (Mean = 0.334ε)

Max = 1.27ε (Mean = 0.334ε)

Max = 1.27ε (Mean = 0.355ε)


Table 8.67. Error rates for ellint_d

GNU C++ version 7.1.0
linux
double

GNU C++ version 7.1.0
linux
long double

Sun compiler version 0x5150
Sun Solaris
long double

Microsoft Visual C++ version 14.1
Win32
double

Elliptic Integral E: Mathworld Data

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 0.862ε (Mean = 0.568ε))

Max = 1.3ε (Mean = 0.813ε)

Max = 1.3ε (Mean = 0.813ε)

Max = 0.862ε (Mean = 0.457ε)

Elliptic Integral D: Random Data

Max = 0ε (Mean = 0ε)

(GSL 2.1: Max = 3.01ε (Mean = 0.928ε))

Max = 2.51ε (Mean = 0.883ε)

Max = 2.51ε (Mean = 0.883ε)

Max = 2.87ε (Mean = 0.805ε)


The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at double precision, and GCC-7.1/Ubuntu for long double and __float128.

Testing

The tests use a mixture of spot test values calculated using values calculated at Wolfram Alpha, and random test data generated using MPFR at 1000-bit precision and a deliberately naive implementation in terms of the Legendre integrals.

Implementation

The implementation for D(φ, k) first performs argument reduction using the relations:

D(-φ, k) = -D(φ, k)

and

D(nπ+φ, k) = 2nD(k) + D(φ, k)

to move φ to the range [0, π/2].

The functions are then implemented in terms of Carlson's integral RD using the relation:


PrevUpHomeNext