]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |