boost/math/special_functions/ellint_rd.hpp

   1 //  Copyright (c) 2006 Xiaogang Zhang
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5 //
   6 //  History:
   7 //  XZ wrote the original of this file as part of the Google
   8 //  Summer of Code 2006.  JM modified it slightly to fit into the
   9 //  Boost.Math conceptual framework better.
  10
  11 #ifndef BOOST_MATH_ELLINT_RD_HPP
  12 #define BOOST_MATH_ELLINT_RD_HPP
  13
  14 #ifdef _MSC_VER
  15 #pragma once
  16 #endif
  17
  18 #include <boost/math/special_functions/math_fwd.hpp>
  19 #include <boost/math/tools/config.hpp>
  20 #include <boost/math/policies/error_handling.hpp>
  21
  22 // Carlson's elliptic integral of the second kind
  23 // R_D(x, y, z) = R_J(x, y, z, z) = 1.5 * \int_{0}^{\infty} [(t+x)(t+y)]^{-1/2} (t+z)^{-3/2} dt
  24 // Carlson, Numerische Mathematik, vol 33, 1 (1979)
  25
  26 namespace boost { namespace math { namespace detail{
  27
  28 template <typename T, typename Policy>
  29 T ellint_rd_imp(T x, T y, T z, const Policy& pol)
  30 {
  31     T value, u, lambda, sigma, factor, tolerance;
  32     T X, Y, Z, EA, EB, EC, ED, EE, S1, S2;
  33     unsigned long k;
  34
  35     BOOST_MATH_STD_USING
  36     using namespace boost::math::tools;
  37
  38     static const char* function = "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)";
  39
  40     if (x < 0)
  41     {
  42        return policies::raise_domain_error<T>(function,
  43             "Argument x must be >= 0, but got %1%", x, pol);
  44     }
  45     if (y < 0)
  46     {
  47        return policies::raise_domain_error<T>(function,
  48             "Argument y must be >= 0, but got %1%", y, pol);
  49     }
  50     if (z <= 0)
  51     {
  52        return policies::raise_domain_error<T>(function,
  53             "Argument z must be > 0, but got %1%", z, pol);
  54     }
  55     if (x + y == 0)
  56     {
  57        return policies::raise_domain_error<T>(function,
  58             "At most one argument can be zero, but got, x + y = %1%", x+y, pol);
  59     }
  60
  61     // error scales as the 6th power of tolerance
  62     tolerance = pow(tools::epsilon<T>() / 3, T(1)/6);
  63
  64     // duplication
  65     sigma = 0;
  66     factor = 1;
  67     k = 1;
  68     do
  69     {
  70         u = (x + y + z + z + z) / 5;
  71         X = (u - x) / u;
  72         Y = (u - y) / u;
  73         Z = (u - z) / u;
  74         if ((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance)
  75            break;
  76         T sx = sqrt(x);
  77         T sy = sqrt(y);
  78         T sz = sqrt(z);
  79         lambda = sy * (sx + sz) + sz * sx; //sqrt(x * y) + sqrt(y * z) + sqrt(z * x);
  80         sigma += factor / (sz * (z + lambda));
  81         factor /= 4;
  82         x = (x + lambda) / 4;
  83         y = (y + lambda) / 4;
  84         z = (z + lambda) / 4;
  85         ++k;
  86     }
  87     while(k < policies::get_max_series_iterations<Policy>());
  88
  89     // Check to see if we gave up too soon:
  90     policies::check_series_iterations<T>(function, k, pol);
  91
  92     // Taylor series expansion to the 5th order
  93     EA = X * Y;
  94     EB = Z * Z;
  95     EC = EA - EB;
  96     ED = EA - 6 * EB;
  97     EE = ED + EC + EC;
  98     S1 = ED * (ED * T(9) / 88 - Z * EE * T(9) / 52 - T(3) / 14);
  99     S2 = Z * (EE / 6 + Z * (-EC * T(9) / 22 + Z * EA * T(3) / 26));
 100     value = 3 * sigma + factor * (1 + S1 + S2) / (u * sqrt(u));
 101
 102     return value;
 103 }
 104
 105 } // namespace detail
 106
 107 template <class T1, class T2, class T3, class Policy>
 108 inline typename tools::promote_args<T1, T2, T3>::type
 109    ellint_rd(T1 x, T2 y, T3 z, const Policy& pol)
 110 {
 111    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
 112    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 113    return policies::checked_narrowing_cast<result_type, Policy>(
 114       detail::ellint_rd_imp(
 115          static_cast<value_type>(x),
 116          static_cast<value_type>(y),
 117          static_cast<value_type>(z), pol), "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)");
 118 }
 119
 120 template <class T1, class T2, class T3>
 121 inline typename tools::promote_args<T1, T2, T3>::type
 122    ellint_rd(T1 x, T2 y, T3 z)
 123 {
 124    return ellint_rd(x, y, z, policies::policy<>());
 125 }
 126
 127 }} // namespaces
 128
 129 #endif // BOOST_MATH_ELLINT_RD_HPP
 130