1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
|
// (C) Copyright John Maddock 2006.
// 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)
#ifndef BOOST_MATH_SF_CBRT_HPP
#define BOOST_MATH_SF_CBRT_HPP
#ifdef _MSC_VER
#pragma once
#endif
#include <boost/math/tools/rational.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/mpl/divides.hpp>
#include <boost/mpl/plus.hpp>
#include <boost/mpl/if.hpp>
#include <boost/type_traits/is_convertible.hpp>
namespace boost{ namespace math{
namespace detail
{
struct big_int_type
{
operator boost::uintmax_t()const;
};
template <class T>
struct largest_cbrt_int_type
{
typedef typename mpl::if_<
boost::is_convertible<big_int_type, T>,
boost::uintmax_t,
unsigned int
>::type type;
};
template <class T, class Policy>
T cbrt_imp(T z, const Policy& pol)
{
BOOST_MATH_STD_USING
//
// cbrt approximation for z in the range [0.5,1]
// It's hard to say what number of terms gives the optimum
// trade off between precision and performance, this seems
// to be about the best for double precision.
//
// Maximum Deviation Found: 1.231e-006
// Expected Error Term: -1.231e-006
// Maximum Relative Change in Control Points: 5.982e-004
//
static const T P[] = {
static_cast<T>(0.37568269008611818),
static_cast<T>(1.3304968705558024),
static_cast<T>(-1.4897101632445036),
static_cast<T>(1.2875573098219835),
static_cast<T>(-0.6398703759826468),
static_cast<T>(0.13584489959258635),
};
static const T correction[] = {
static_cast<T>(0.62996052494743658238360530363911), // 2^-2/3
static_cast<T>(0.79370052598409973737585281963615), // 2^-1/3
static_cast<T>(1),
static_cast<T>(1.2599210498948731647672106072782), // 2^1/3
static_cast<T>(1.5874010519681994747517056392723), // 2^2/3
};
if((boost::math::isinf)(z) || (z == 0))
return z;
if(!(boost::math::isfinite)(z))
{
return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol);
}
int i_exp, sign(1);
if(z < 0)
{
z = -z;
sign = -sign;
}
T guess = frexp(z, &i_exp);
int original_i_exp = i_exp; // save for later
guess = tools::evaluate_polynomial(P, guess);
int i_exp3 = i_exp / 3;
typedef typename largest_cbrt_int_type<T>::type shift_type;
BOOST_STATIC_ASSERT( ::std::numeric_limits<shift_type>::radix == 2);
if(abs(i_exp3) < std::numeric_limits<shift_type>::digits)
{
if(i_exp3 > 0)
guess *= shift_type(1u) << i_exp3;
else
guess /= shift_type(1u) << -i_exp3;
}
else
{
guess = ldexp(guess, i_exp3);
}
i_exp %= 3;
guess *= correction[i_exp + 2];
//
// Now inline Halley iteration.
// We do this here rather than calling tools::halley_iterate since we can
// simplify the expressions algebraically, and don't need most of the error
// checking of the boilerplate version as we know in advance that the function
// is well behaved...
//
typedef typename policies::precision<T, Policy>::type prec;
typedef typename mpl::divides<prec, mpl::int_<3> >::type prec3;
typedef typename mpl::plus<prec3, mpl::int_<3> >::type new_prec;
typedef typename policies::normalise<Policy, policies::digits2<new_prec::value> >::type new_policy;
//
// Epsilon calculation uses compile time arithmetic when it's available for type T,
// otherwise uses ldexp to calculate at runtime:
//
T eps = (new_prec::value > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3);
T diff;
if(original_i_exp < std::numeric_limits<T>::max_exponent - 3)
{
//
// Safe from overflow, use the fast method:
//
do
{
T g3 = guess * guess * guess;
diff = (g3 + z + z) / (g3 + g3 + z);
guess *= diff;
}
while(fabs(1 - diff) > eps);
}
else
{
//
// Either we're ready to overflow, or we can't tell because numeric_limits isn't
// available for type T:
//
do
{
T g2 = guess * guess;
diff = (g2 - z / guess) / (2 * guess + z / g2);
guess -= diff;
}
while((guess * eps) < fabs(diff));
}
return sign * guess;
}
} // namespace detail
template <class T, class Policy>
inline typename tools::promote_args<T>::type cbrt(T z, const Policy& pol)
{
typedef typename tools::promote_args<T>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
return static_cast<result_type>(detail::cbrt_imp(value_type(z), pol));
}
template <class T>
inline typename tools::promote_args<T>::type cbrt(T z)
{
return cbrt(z, policies::policy<>());
}
} // namespace math
} // namespace boost
#endif // BOOST_MATH_SF_CBRT_HPP
|
