ellint_rf.hpp 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168
  1. // Copyright (c) 2006 Xiaogang Zhang, 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. // History:
  7. // XZ wrote the original of this file as part of the Google
  8. // Summer of Code 2006. JM modified it to fit into the
  9. // Boost.Math conceptual framework better, and to handle
  10. // types longer than 80-bit reals.
  11. // Updated 2015 to use Carlson's latest methods.
  12. //
  13. #ifndef BOOST_MATH_ELLINT_RF_HPP
  14. #define BOOST_MATH_ELLINT_RF_HPP
  15. #ifdef _MSC_VER
  16. #pragma once
  17. #endif
  18. #include <boost/math/special_functions/math_fwd.hpp>
  19. #include <boost/math/tools/config.hpp>
  20. #include <boost/math/constants/constants.hpp>
  21. #include <boost/math/policies/error_handling.hpp>
  22. #include <boost/math/special_functions/ellint_rc.hpp>
  23. // Carlson's elliptic integral of the first kind
  24. // R_F(x, y, z) = 0.5 * \int_{0}^{\infty} [(t+x)(t+y)(t+z)]^{-1/2} dt
  25. // Carlson, Numerische Mathematik, vol 33, 1 (1979)
  26. namespace boost { namespace math { namespace detail{
  27. template <typename T, typename Policy>
  28. T ellint_rf_imp(T x, T y, T z, const Policy& pol)
  29. {
  30. BOOST_MATH_STD_USING
  31. using namespace boost::math;
  32. using std::swap;
  33. static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)";
  34. if(x < 0 || y < 0 || z < 0)
  35. {
  36. return policies::raise_domain_error<T>(function, "domain error, all arguments must be non-negative, only sensible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
  37. }
  38. if(x + y == 0 || y + z == 0 || z + x == 0)
  39. {
  40. return policies::raise_domain_error<T>(function, "domain error, at most one argument can be zero, only sensible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
  41. }
  42. //
  43. // Special cases from http://dlmf.nist.gov/19.20#i
  44. //
  45. if(x == y)
  46. {
  47. if(x == z)
  48. {
  49. // x, y, z equal:
  50. return 1 / sqrt(x);
  51. }
  52. else
  53. {
  54. // 2 equal, x and y:
  55. if(z == 0)
  56. return constants::pi<T>() / (2 * sqrt(x));
  57. else
  58. return ellint_rc_imp(z, x, pol);
  59. }
  60. }
  61. if(x == z)
  62. {
  63. if(y == 0)
  64. return constants::pi<T>() / (2 * sqrt(x));
  65. else
  66. return ellint_rc_imp(y, x, pol);
  67. }
  68. if(y == z)
  69. {
  70. if(x == 0)
  71. return constants::pi<T>() / (2 * sqrt(y));
  72. else
  73. return ellint_rc_imp(x, y, pol);
  74. }
  75. if(x == 0)
  76. swap(x, z);
  77. else if(y == 0)
  78. swap(y, z);
  79. if(z == 0)
  80. {
  81. //
  82. // Special case for one value zero:
  83. //
  84. T xn = sqrt(x);
  85. T yn = sqrt(y);
  86. while(fabs(xn - yn) >= T(2.7) * tools::root_epsilon<T>() * fabs(xn))
  87. {
  88. T t = sqrt(xn * yn);
  89. xn = (xn + yn) / 2;
  90. yn = t;
  91. }
  92. return constants::pi<T>() / (xn + yn);
  93. }
  94. T xn = x;
  95. T yn = y;
  96. T zn = z;
  97. T An = (x + y + z) / 3;
  98. T A0 = An;
  99. T Q = pow(3 * boost::math::tools::epsilon<T>(), T(-1) / 8) * (std::max)((std::max)(fabs(An - xn), fabs(An - yn)), fabs(An - zn));
  100. T fn = 1;
  101. // duplication
  102. unsigned k = 1;
  103. for(; k < boost::math::policies::get_max_series_iterations<Policy>(); ++k)
  104. {
  105. T root_x = sqrt(xn);
  106. T root_y = sqrt(yn);
  107. T root_z = sqrt(zn);
  108. T lambda = root_x * root_y + root_x * root_z + root_y * root_z;
  109. An = (An + lambda) / 4;
  110. xn = (xn + lambda) / 4;
  111. yn = (yn + lambda) / 4;
  112. zn = (zn + lambda) / 4;
  113. Q /= 4;
  114. fn *= 4;
  115. if(Q < fabs(An))
  116. break;
  117. }
  118. // Check to see if we gave up too soon:
  119. policies::check_series_iterations<T>(function, k, pol);
  120. BOOST_MATH_INSTRUMENT_VARIABLE(k);
  121. T X = (A0 - x) / (An * fn);
  122. T Y = (A0 - y) / (An * fn);
  123. T Z = -X - Y;
  124. // Taylor series expansion to the 7th order
  125. T E2 = X * Y - Z * Z;
  126. T E3 = X * Y * Z;
  127. return (1 + E3 * (T(1) / 14 + 3 * E3 / 104) + E2 * (T(-1) / 10 + E2 / 24 - (3 * E3) / 44 - 5 * E2 * E2 / 208 + E2 * E3 / 16)) / sqrt(An);
  128. }
  129. } // namespace detail
  130. template <class T1, class T2, class T3, class Policy>
  131. inline typename tools::promote_args<T1, T2, T3>::type
  132. ellint_rf(T1 x, T2 y, T3 z, const Policy& pol)
  133. {
  134. typedef typename tools::promote_args<T1, T2, T3>::type result_type;
  135. typedef typename policies::evaluation<result_type, Policy>::type value_type;
  136. return policies::checked_narrowing_cast<result_type, Policy>(
  137. detail::ellint_rf_imp(
  138. static_cast<value_type>(x),
  139. static_cast<value_type>(y),
  140. static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)");
  141. }
  142. template <class T1, class T2, class T3>
  143. inline typename tools::promote_args<T1, T2, T3>::type
  144. ellint_rf(T1 x, T2 y, T3 z)
  145. {
  146. return ellint_rf(x, y, z, policies::policy<>());
  147. }
  148. }} // namespaces
  149. #endif // BOOST_MATH_ELLINT_RF_HPP