|
Home | Libraries | People | FAQ | More |
#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
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.
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.
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 |
GNU C++ version 5.1.0 |
GNU C++ version 5.1.0 |
Sun compiler version 0x5130 |
|
|---|---|---|---|---|
|
RC: Random data |
Max = 0.962ε (Mean = 0.407ε) |
Max = 0ε (Mean = 0ε) |
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 |
GNU C++ version 5.1.0 |
GNU C++ version 5.1.0 |
Sun compiler version 0x5130 |
|
|---|---|---|---|---|
|
RD: Random data |
Max = 2.16ε (Mean = 0.803ε) |
Max = 0ε (Mean = 0ε) |
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ε) |
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ε) |
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ε) |
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ε) |
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ε) |
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 |
GNU C++ version 5.1.0 |
GNU C++ version 5.1.0 |
Sun compiler version 0x5130 |
|
|---|---|---|---|---|
|
RG: Random Data |
Max = 3.65ε (Mean = 0.929ε) |
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ε) |
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ε) |
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ε) |
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ε) |
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 |
GNU C++ version 5.1.0 |
GNU C++ version 5.1.0 |
Sun compiler version 0x5130 |
|
|---|---|---|---|---|
|
RF: Random data |
Max = 2.02ε (Mean = 0.677ε) |
Max = 0ε (Mean = 0ε) |
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ε) |
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ε) |
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ε) |
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ε) |
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 |
GNU C++ version 5.1.0 |
GNU C++ version 5.1.0 |
Sun compiler version 0x5130 |
|
|---|---|---|---|---|
|
RJ: Random data |
Max = 119ε (Mean = 4.32ε) |
Max = 0.52ε (Mean = 0.0184ε) |
Max = 186ε (Mean = 6.67ε) |
Max = 186ε (Mean = 6.7ε) |
|
RJ: 4 Equal Values |
Max = 1.03ε (Mean = 0.418ε) |
Max = 0ε (Mean = 0ε) |
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ε) |
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ε) |
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ε) |
Max = 17.2ε (Mean = 1.16ε) |
Max = 16.6ε (Mean = 1.15ε) |
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: