| Home | Libraries | People | FAQ | More

### 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: