ceph/src/boost/libs/math/include/boost/math/special_functions/ellint_rg.hpp

   1 //  Copyright (c) 2015 John Maddock
   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 #ifndef BOOST_MATH_ELLINT_RG_HPP
   7 #define BOOST_MATH_ELLINT_RG_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/math_fwd.hpp>
  14 #include <boost/math/tools/config.hpp>
  15 #include <boost/math/constants/constants.hpp>
  16 #include <boost/math/policies/error_handling.hpp>
  17 #include <boost/math/special_functions/ellint_rd.hpp>
  18 #include <boost/math/special_functions/ellint_rf.hpp>
  19 #include <boost/math/special_functions/pow.hpp>
  20
  21 namespace boost { namespace math { namespace detail{
  22
  23    template <typename T, typename Policy>
  24    T ellint_rg_imp(T x, T y, T z, const Policy& pol)
  25    {
  26       BOOST_MATH_STD_USING
  27       static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)";
  28
  29       if(x < 0 || y < 0 || z < 0)
  30       {
  31          return policies::raise_domain_error<T>(function,
  32             "domain error, all arguments must be non-negative, "
  33             "only sensible result is %1%.",
  34             std::numeric_limits<T>::quiet_NaN(), pol);
  35       }
  36       //
  37       // Function is symmetric in x, y and z, but we require
  38       // (x - z)(y - z) >= 0 to avoid cancellation error in the result
  39       // which implies (for example) x >= z >= y
  40       //
  41       using std::swap;
  42       if(x < y)
  43          swap(x, y);
  44       if(x < z)
  45          swap(x, z);
  46       if(y > z)
  47          swap(y, z);
  48
  49       BOOST_ASSERT(x >= z);
  50       BOOST_ASSERT(z >= y);
  51       //
  52       // Special cases from http://dlmf.nist.gov/19.20#ii
  53       //
  54       if(x == z)
  55       {
  56          if(y == z)
  57          {
  58             // x = y = z
  59             // This also works for x = y = z = 0 presumably.
  60             return sqrt(x);
  61          }
  62          else if(y == 0)
  63          {
  64             // x = y, z = 0
  65             return constants::pi<T>() * sqrt(x) / 4;
  66          }
  67          else
  68          {
  69             // x = z, y != 0
  70             swap(x, y);
  71             return (x == 0) ? T(sqrt(z) / 2) : T((z * ellint_rc_imp(x, z, pol) + sqrt(x)) / 2);
  72          }
  73       }
  74       else if(y == z)
  75       {
  76          if(x == 0)
  77             return constants::pi<T>() * sqrt(y) / 4;
  78          else
  79             return (y == 0) ? T(sqrt(x) / 2) : T((y * ellint_rc_imp(x, y, pol) + sqrt(x)) / 2);
  80       }
  81       else if(y == 0)
  82       {
  83          swap(y, z);
  84          //
  85          // Special handling for common case, from
  86          // Numerical Computation of Real or Complex Elliptic Integrals, eq.46
  87          //
  88          T xn = sqrt(x);
  89          T yn = sqrt(y);
  90          T x0 = xn;
  91          T y0 = yn;
  92          T sum = 0;
  93          T sum_pow = 0.25f;
  94
  95          while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn))
  96          {
  97             T t = sqrt(xn * yn);
  98             xn = (xn + yn) / 2;
  99             yn = t;
 100             sum_pow *= 2;
 101             sum += sum_pow * boost::math::pow<2>(xn - yn);
 102          }
 103          T RF = constants::pi<T>() / (xn + yn);
 104          return ((boost::math::pow<2>((x0 + y0) / 2) - sum) * RF) / 2;
 105       }
 106       return (z * ellint_rf_imp(x, y, z, pol)
 107          - (x - z) * (y - z) * ellint_rd_imp(x, y, z, pol) / 3
 108          + sqrt(x * y / z)) / 2;
 109    }
 110
 111 } // namespace detail
 112
 113 template <class T1, class T2, class T3, class Policy>
 114 inline typename tools::promote_args<T1, T2, T3>::type
 115    ellint_rg(T1 x, T2 y, T3 z, const Policy& pol)
 116 {
 117    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
 118    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 119    return policies::checked_narrowing_cast<result_type, Policy>(
 120       detail::ellint_rg_imp(
 121          static_cast<value_type>(x),
 122          static_cast<value_type>(y),
 123          static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)");
 124 }
 125
 126 template <class T1, class T2, class T3>
 127 inline typename tools::promote_args<T1, T2, T3>::type
 128    ellint_rg(T1 x, T2 y, T3 z)
 129 {
 130    return ellint_rg(x, y, z, policies::policy<>());
 131 }
 132
 133 }} // namespaces
 134
 135 #endif // BOOST_MATH_ELLINT_RG_HPP
 136