ellint_carlson.qbk

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

[/
Copyright (c) 2006 Xiaogang Zhang
Copyright (c) 2006 John Maddock
Use, modification and distribution are subject to the
Boost Software License, Version 1.0. (See accompanying file
LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
]

[section:ellint_carlson Elliptic Integrals - Carlson Form]

[heading Synopsis]

``
  #include <boost/math/special_functions/ellint_rf.hpp>
``

  namespace boost { namespace math {

  template <class T1, class T2, class T3>
  ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z)

  template <class T1, class T2, class T3, class ``__Policy``>
  ``__sf_result`` 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>
  ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z)

  template <class T1, class T2, class T3, class ``__Policy``>
  ``__sf_result`` 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>
  ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p)

  template <class T1, class T2, class T3, class T4, class ``__Policy``>
  ``__sf_result`` 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>
  ``__sf_result`` ellint_rc(T1 x, T2 y)

  template <class T1, class T2, class ``__Policy``>
  ``__sf_result`` ellint_rc(T1 x, T2 y, const ``__Policy``&)

  }} // namespaces


[heading 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:

[graph ellint_carlson]

The return type of these functions is computed using the __arg_pomotion_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>
  ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z)
  
  template <class T1, class T2, class T3, class ``__Policy``>
  ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z, const ``__Policy``&)
  
Returns Carlson's Elliptic Integral R[sub F]:

[equation ellint9]

Requires that all of the arguments are non-negative, and at most
one may be zero.  Otherwise returns the result of __domain_error.

[optional_policy]

  template <class T1, class T2, class T3>
  ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z)
  
  template <class T1, class T2, class T3, class ``__Policy``>
  ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z, const ``__Policy``&)
  
Returns Carlson's elliptic integral R[sub D]:

[equation ellint10]

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.

[optional_policy]

  template <class T1, class T2, class T3, class T4>
  ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p)

  template <class T1, class T2, class T3, class T4, class ``__Policy``>
  ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p, const ``__Policy``&)

Returns Carlson's elliptic integral R[sub J]:
  
[equation ellint11]

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.

[optional_policy]

When ['p < 0] the function returns the
[@http://en.wikipedia.org/wiki/Cauchy_principal_value Cauchy principal value]
using the relation:

[equation ellint17]

  template <class T1, class T2>
  ``__sf_result`` ellint_rc(T1 x, T2 y)

  template <class T1, class T2, class ``__Policy``>
  ``__sf_result`` ellint_rc(T1 x, T2 y, const ``__Policy``&)

Returns Carlson's elliptic integral R[sub C]:
  
[equation ellint12]

Requires that ['x > 0] and that ['y != 0].  
Otherwise returns the result of __domain_error.

[optional_policy]

When ['y < 0] the function returns the
[@http://mathworld.wolfram.com/CauchyPrincipalValue.html Cauchy principal value]
using the relation:

[equation ellint18]

[heading Testing]

There are two sets of tests.

Spot tests compare selected values with test data given in:

[:B. C. Carlson, ['[@http://arxiv.org/abs/math.CA/9409227 
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.

[heading 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 __zero_error.  All values
are relative errors in units of epsilon.

[table Errors Rates in the Carlson Elliptic Integrals
[[Significand Size] [Platform and Compiler] [R[sub F]] [R[sub D]] [R[sub J]] [R[sub C]]]
[[53] [Win32 / Visual C++ 8.0] [Peak=2.9 Mean=0.75] [Peak=2.6 Mean=0.9] [Peak=108 Mean=6.9] [Peak=2.4 Mean=0.6] ]
[[64] [Red Hat Linux / G++ 3.4] [Peak=2.5 Mean=0.75] [Peak=2.7 Mean=0.9] [Peak=105 Mean=8] [Peak=1.9 Mean=0.7] ]
[[113] [HP-UX / HP aCC 6] [Peak=5.3 Mean=1.6] [Peak=2.9 Mean=0.99] [Peak=180 Mean=12] [Peak=1.8 Mean=0.7] ]
]

[heading Implementation]

The key of Carlson's algorithm [[link ellint_ref_carlson79 Carlson79]] is the
duplication theorem:

[equation ellint13]

By applying it repeatedly, ['x], ['y], ['z] get
closer and closer. When they are nearly equal, the special case equation

[equation ellint16]

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.

For ['p < 0] in ['R[sub J](x, y, z, p)] and ['y < 0] in ['R[sub C](x, y)],
the integrals are singular and their
[@http://mathworld.wolfram.com/CauchyPrincipalValue.html Cauchy principal values]
are returned via the relations:

[equation ellint17]

[equation ellint18]

[endsect]