]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | |
2 | // Copyright Christopher Kormanyos 2002 - 2013. | |
3 | // Copyright 2011 - 2013 John Maddock. Distributed under the Boost | |
4 | // Distributed under the Boost Software License, Version 1.0. | |
5 | // (See accompanying file LICENSE_1_0.txt or copy at | |
6 | // http://www.boost.org/LICENSE_1_0.txt) | |
7 | ||
8 | // This work is based on an earlier work: | |
9 | // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations", | |
10 | // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469 | |
11 | // | |
12 | // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp | |
92f5a8d4 | 13 | // |
7c673cae FG |
14 | |
15 | #ifdef BOOST_MSVC | |
16 | #pragma warning(push) | |
92f5a8d4 | 17 | #pragma warning(disable : 6326) // comparison of two constants |
7c673cae FG |
18 | #endif |
19 | ||
92f5a8d4 | 20 | namespace detail { |
7c673cae | 21 | |
92f5a8d4 | 22 | template <typename T, typename U> |
7c673cae FG |
23 | inline void pow_imp(T& result, const T& t, const U& p, const mpl::false_&) |
24 | { | |
25 | // Compute the pure power of typename T t^p. | |
26 | // Use the S-and-X binary method, as described in | |
27 | // D. E. Knuth, "The Art of Computer Programming", Vol. 2, | |
28 | // Section 4.6.3 . The resulting computational complexity | |
29 | // is order log2[abs(p)]. | |
30 | ||
31 | typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type; | |
32 | ||
92f5a8d4 | 33 | if (&result == &t) |
7c673cae FG |
34 | { |
35 | T temp; | |
36 | pow_imp(temp, t, p, mpl::false_()); | |
37 | result = temp; | |
38 | return; | |
39 | } | |
40 | ||
41 | // This will store the result. | |
92f5a8d4 | 42 | if (U(p % U(2)) != U(0)) |
7c673cae FG |
43 | { |
44 | result = t; | |
45 | } | |
46 | else | |
47 | result = int_type(1); | |
48 | ||
49 | U p2(p); | |
50 | ||
51 | // The variable x stores the binary powers of t. | |
52 | T x(t); | |
53 | ||
92f5a8d4 | 54 | while (U(p2 /= 2) != U(0)) |
7c673cae FG |
55 | { |
56 | // Square x for each binary power. | |
57 | eval_multiply(x, x); | |
58 | ||
59 | const bool has_binary_power = (U(p2 % U(2)) != U(0)); | |
60 | ||
92f5a8d4 | 61 | if (has_binary_power) |
7c673cae FG |
62 | { |
63 | // Multiply the result with each binary power contained in the exponent. | |
64 | eval_multiply(result, x); | |
65 | } | |
66 | } | |
67 | } | |
68 | ||
92f5a8d4 | 69 | template <typename T, typename U> |
7c673cae FG |
70 | inline void pow_imp(T& result, const T& t, const U& p, const mpl::true_&) |
71 | { | |
72 | // Signed integer power, just take care of the sign then call the unsigned version: | |
92f5a8d4 TL |
73 | typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type; |
74 | typedef typename make_unsigned<U>::type ui_type; | |
7c673cae | 75 | |
92f5a8d4 | 76 | if (p < 0) |
7c673cae FG |
77 | { |
78 | T temp; | |
79 | temp = static_cast<int_type>(1); | |
80 | T denom; | |
81 | pow_imp(denom, t, static_cast<ui_type>(-p), mpl::false_()); | |
82 | eval_divide(result, temp, denom); | |
83 | return; | |
84 | } | |
85 | pow_imp(result, t, static_cast<ui_type>(p), mpl::false_()); | |
86 | } | |
87 | ||
88 | } // namespace detail | |
89 | ||
92f5a8d4 TL |
90 | template <typename T, typename U> |
91 | inline typename enable_if_c<is_integral<U>::value>::type eval_pow(T& result, const T& t, const U& p) | |
7c673cae FG |
92 | { |
93 | detail::pow_imp(result, t, p, boost::is_signed<U>()); | |
94 | } | |
95 | ||
96 | template <class T> | |
97 | void hyp0F0(T& H0F0, const T& x) | |
98 | { | |
99 | // Compute the series representation of Hypergeometric0F0 taken from | |
100 | // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F0/06/01/ | |
101 | // There are no checks on input range or parameter boundaries. | |
102 | ||
103 | typedef typename mpl::front<typename T::unsigned_types>::type ui_type; | |
104 | ||
105 | BOOST_ASSERT(&H0F0 != &x); | |
106 | long tol = boost::multiprecision::detail::digits2<number<T, et_on> >::value(); | |
92f5a8d4 | 107 | T t; |
7c673cae FG |
108 | |
109 | T x_pow_n_div_n_fact(x); | |
110 | ||
111 | eval_add(H0F0, x_pow_n_div_n_fact, ui_type(1)); | |
112 | ||
113 | T lim; | |
114 | eval_ldexp(lim, H0F0, 1 - tol); | |
92f5a8d4 | 115 | if (eval_get_sign(lim) < 0) |
7c673cae FG |
116 | lim.negate(); |
117 | ||
118 | ui_type n; | |
119 | ||
92f5a8d4 TL |
120 | const unsigned series_limit = |
121 | boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100 | |
122 | ? 100 | |
123 | : boost::multiprecision::detail::digits2<number<T, et_on> >::value(); | |
7c673cae | 124 | // Series expansion of hyperg_0f0(; ; x). |
92f5a8d4 | 125 | for (n = 2; n < series_limit; ++n) |
7c673cae FG |
126 | { |
127 | eval_multiply(x_pow_n_div_n_fact, x); | |
128 | eval_divide(x_pow_n_div_n_fact, n); | |
129 | eval_add(H0F0, x_pow_n_div_n_fact); | |
130 | bool neg = eval_get_sign(x_pow_n_div_n_fact) < 0; | |
92f5a8d4 | 131 | if (neg) |
7c673cae | 132 | x_pow_n_div_n_fact.negate(); |
92f5a8d4 | 133 | if (lim.compare(x_pow_n_div_n_fact) > 0) |
7c673cae | 134 | break; |
92f5a8d4 | 135 | if (neg) |
7c673cae FG |
136 | x_pow_n_div_n_fact.negate(); |
137 | } | |
92f5a8d4 | 138 | if (n >= series_limit) |
7c673cae FG |
139 | BOOST_THROW_EXCEPTION(std::runtime_error("H0F0 failed to converge")); |
140 | } | |
141 | ||
142 | template <class T> | |
143 | void hyp1F0(T& H1F0, const T& a, const T& x) | |
144 | { | |
145 | // Compute the series representation of Hypergeometric1F0 taken from | |
146 | // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F0/06/01/01/ | |
147 | // and also see the corresponding section for the power function (i.e. x^a). | |
148 | // There are no checks on input range or parameter boundaries. | |
149 | ||
150 | typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type; | |
151 | ||
152 | BOOST_ASSERT(&H1F0 != &x); | |
153 | BOOST_ASSERT(&H1F0 != &a); | |
154 | ||
155 | T x_pow_n_div_n_fact(x); | |
92f5a8d4 TL |
156 | T pochham_a(a); |
157 | T ap(a); | |
7c673cae FG |
158 | |
159 | eval_multiply(H1F0, pochham_a, x_pow_n_div_n_fact); | |
160 | eval_add(H1F0, si_type(1)); | |
161 | T lim; | |
162 | eval_ldexp(lim, H1F0, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value()); | |
92f5a8d4 | 163 | if (eval_get_sign(lim) < 0) |
7c673cae FG |
164 | lim.negate(); |
165 | ||
166 | si_type n; | |
92f5a8d4 | 167 | T term, part; |
7c673cae FG |
168 | |
169 | const si_type series_limit = | |
92f5a8d4 TL |
170 | boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100 |
171 | ? 100 | |
172 | : boost::multiprecision::detail::digits2<number<T, et_on> >::value(); | |
7c673cae | 173 | // Series expansion of hyperg_1f0(a; ; x). |
92f5a8d4 | 174 | for (n = 2; n < series_limit; n++) |
7c673cae FG |
175 | { |
176 | eval_multiply(x_pow_n_div_n_fact, x); | |
177 | eval_divide(x_pow_n_div_n_fact, n); | |
178 | eval_increment(ap); | |
179 | eval_multiply(pochham_a, ap); | |
180 | eval_multiply(term, pochham_a, x_pow_n_div_n_fact); | |
181 | eval_add(H1F0, term); | |
92f5a8d4 | 182 | if (eval_get_sign(term) < 0) |
7c673cae | 183 | term.negate(); |
92f5a8d4 | 184 | if (lim.compare(term) >= 0) |
7c673cae FG |
185 | break; |
186 | } | |
92f5a8d4 | 187 | if (n >= series_limit) |
7c673cae FG |
188 | BOOST_THROW_EXCEPTION(std::runtime_error("H1F0 failed to converge")); |
189 | } | |
190 | ||
191 | template <class T> | |
192 | void eval_exp(T& result, const T& x) | |
193 | { | |
194 | BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The exp function is only valid for floating point types."); | |
92f5a8d4 | 195 | if (&x == &result) |
7c673cae FG |
196 | { |
197 | T temp; | |
198 | eval_exp(temp, x); | |
199 | result = temp; | |
200 | return; | |
201 | } | |
202 | typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; | |
92f5a8d4 TL |
203 | typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type; |
204 | typedef typename T::exponent_type exp_type; | |
7c673cae FG |
205 | typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type; |
206 | ||
207 | // Handle special arguments. | |
92f5a8d4 | 208 | int type = eval_fpclassify(x); |
7c673cae | 209 | bool isneg = eval_get_sign(x) < 0; |
92f5a8d4 | 210 | if (type == (int)FP_NAN) |
7c673cae FG |
211 | { |
212 | result = x; | |
92f5a8d4 | 213 | errno = EDOM; |
7c673cae FG |
214 | return; |
215 | } | |
92f5a8d4 | 216 | else if (type == (int)FP_INFINITE) |
7c673cae | 217 | { |
92f5a8d4 | 218 | if (isneg) |
7c673cae | 219 | result = ui_type(0u); |
92f5a8d4 | 220 | else |
7c673cae FG |
221 | result = x; |
222 | return; | |
223 | } | |
92f5a8d4 | 224 | else if (type == (int)FP_ZERO) |
7c673cae FG |
225 | { |
226 | result = ui_type(1); | |
227 | return; | |
228 | } | |
229 | ||
230 | // Get local copy of argument and force it to be positive. | |
231 | T xx = x; | |
232 | T exp_series; | |
92f5a8d4 | 233 | if (isneg) |
7c673cae FG |
234 | xx.negate(); |
235 | ||
236 | // Check the range of the argument. | |
92f5a8d4 | 237 | if (xx.compare(si_type(1)) <= 0) |
7c673cae FG |
238 | { |
239 | // | |
240 | // Use series for exp(x) - 1: | |
241 | // | |
242 | T lim; | |
92f5a8d4 | 243 | if (std::numeric_limits<number<T, et_on> >::is_specialized) |
7c673cae FG |
244 | lim = std::numeric_limits<number<T, et_on> >::epsilon().backend(); |
245 | else | |
246 | { | |
247 | result = ui_type(1); | |
248 | eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value()); | |
249 | } | |
250 | unsigned k = 2; | |
251 | exp_series = xx; | |
92f5a8d4 TL |
252 | result = si_type(1); |
253 | if (isneg) | |
7c673cae FG |
254 | eval_subtract(result, exp_series); |
255 | else | |
256 | eval_add(result, exp_series); | |
257 | eval_multiply(exp_series, xx); | |
258 | eval_divide(exp_series, ui_type(k)); | |
259 | eval_add(result, exp_series); | |
92f5a8d4 | 260 | while (exp_series.compare(lim) > 0) |
7c673cae FG |
261 | { |
262 | ++k; | |
263 | eval_multiply(exp_series, xx); | |
264 | eval_divide(exp_series, ui_type(k)); | |
92f5a8d4 | 265 | if (isneg && (k & 1)) |
7c673cae FG |
266 | eval_subtract(result, exp_series); |
267 | else | |
268 | eval_add(result, exp_series); | |
269 | } | |
270 | return; | |
271 | } | |
272 | ||
273 | // Check for pure-integer arguments which can be either signed or unsigned. | |
274 | typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type ll; | |
275 | eval_trunc(exp_series, x); | |
276 | eval_convert_to(&ll, exp_series); | |
92f5a8d4 | 277 | if (x.compare(ll) == 0) |
7c673cae FG |
278 | { |
279 | detail::pow_imp(result, get_constant_e<T>(), ll, mpl::true_()); | |
280 | return; | |
281 | } | |
92f5a8d4 | 282 | else if (exp_series.compare(x) == 0) |
b32b8144 | 283 | { |
92f5a8d4 | 284 | // We have a value that has no fractional part, but is too large to fit |
b32b8144 FG |
285 | // in a long long, in this situation the code below will fail, so |
286 | // we're just going to assume that this will overflow: | |
92f5a8d4 | 287 | if (isneg) |
b32b8144 FG |
288 | result = ui_type(0); |
289 | else | |
290 | result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend(); | |
291 | return; | |
292 | } | |
7c673cae FG |
293 | |
294 | // The algorithm for exp has been taken from MPFUN. | |
295 | // exp(t) = [ (1 + r + r^2/2! + r^3/3! + r^4/4! ...)^p2 ] * 2^n | |
296 | // where p2 is a power of 2 such as 2048, r = t_prime / p2, and | |
297 | // t_prime = t - n*ln2, with n chosen to minimize the absolute | |
298 | // value of t_prime. In the resulting Taylor series, which is | |
299 | // implemented as a hypergeometric function, |r| is bounded by | |
300 | // ln2 / p2. For small arguments, no scaling is done. | |
301 | ||
302 | // Compute the exponential series of the (possibly) scaled argument. | |
303 | ||
304 | eval_divide(result, xx, get_constant_ln2<T>()); | |
305 | exp_type n; | |
306 | eval_convert_to(&n, result); | |
307 | ||
b32b8144 FG |
308 | if (n == (std::numeric_limits<exp_type>::max)()) |
309 | { | |
310 | // Exponent is too large to fit in our exponent type: | |
311 | if (isneg) | |
312 | result = ui_type(0); | |
313 | else | |
314 | result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend(); | |
315 | return; | |
316 | } | |
317 | ||
7c673cae FG |
318 | // The scaling is 2^11 = 2048. |
319 | const si_type p2 = static_cast<si_type>(si_type(1) << 11); | |
320 | ||
321 | eval_multiply(exp_series, get_constant_ln2<T>(), static_cast<canonical_exp_type>(n)); | |
322 | eval_subtract(exp_series, xx); | |
323 | eval_divide(exp_series, p2); | |
324 | exp_series.negate(); | |
325 | hyp0F0(result, exp_series); | |
326 | ||
327 | detail::pow_imp(exp_series, result, p2, mpl::true_()); | |
328 | result = ui_type(1); | |
329 | eval_ldexp(result, result, n); | |
330 | eval_multiply(exp_series, result); | |
331 | ||
92f5a8d4 | 332 | if (isneg) |
7c673cae FG |
333 | eval_divide(result, ui_type(1), exp_series); |
334 | else | |
335 | result = exp_series; | |
336 | } | |
337 | ||
338 | template <class T> | |
339 | void eval_log(T& result, const T& arg) | |
340 | { | |
341 | BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types."); | |
342 | // | |
343 | // We use a variation of http://dlmf.nist.gov/4.45#i | |
344 | // using frexp to reduce the argument to x * 2^n, | |
345 | // then let y = x - 1 and compute: | |
346 | // log(x) = log(2) * n + log1p(1 + y) | |
347 | // | |
348 | typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; | |
92f5a8d4 | 349 | typedef typename T::exponent_type exp_type; |
7c673cae | 350 | typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type; |
92f5a8d4 TL |
351 | typedef typename mpl::front<typename T::float_types>::type fp_type; |
352 | int s = eval_signbit(arg); | |
353 | switch (eval_fpclassify(arg)) | |
b32b8144 FG |
354 | { |
355 | case FP_NAN: | |
356 | result = arg; | |
92f5a8d4 | 357 | errno = EDOM; |
b32b8144 FG |
358 | return; |
359 | case FP_INFINITE: | |
92f5a8d4 TL |
360 | if (s) |
361 | break; | |
b32b8144 FG |
362 | result = arg; |
363 | return; | |
364 | case FP_ZERO: | |
365 | result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend(); | |
366 | result.negate(); | |
367 | errno = ERANGE; | |
368 | return; | |
369 | } | |
92f5a8d4 | 370 | if (s) |
b32b8144 FG |
371 | { |
372 | result = std::numeric_limits<number<T> >::quiet_NaN().backend(); | |
92f5a8d4 | 373 | errno = EDOM; |
b32b8144 FG |
374 | return; |
375 | } | |
7c673cae FG |
376 | |
377 | exp_type e; | |
92f5a8d4 | 378 | T t; |
7c673cae FG |
379 | eval_frexp(t, arg, &e); |
380 | bool alternate = false; | |
381 | ||
92f5a8d4 | 382 | if (t.compare(fp_type(2) / fp_type(3)) <= 0) |
7c673cae FG |
383 | { |
384 | alternate = true; | |
385 | eval_ldexp(t, t, 1); | |
386 | --e; | |
387 | } | |
92f5a8d4 | 388 | |
7c673cae FG |
389 | eval_multiply(result, get_constant_ln2<T>(), canonical_exp_type(e)); |
390 | INSTRUMENT_BACKEND(result); | |
391 | eval_subtract(t, ui_type(1)); /* -0.3 <= t <= 0.3 */ | |
92f5a8d4 | 392 | if (!alternate) |
7c673cae FG |
393 | t.negate(); /* 0 <= t <= 0.33333 */ |
394 | T pow = t; | |
395 | T lim; | |
396 | T t2; | |
397 | ||
92f5a8d4 | 398 | if (alternate) |
7c673cae FG |
399 | eval_add(result, t); |
400 | else | |
401 | eval_subtract(result, t); | |
402 | ||
92f5a8d4 | 403 | if (std::numeric_limits<number<T, et_on> >::is_specialized) |
7c673cae FG |
404 | eval_multiply(lim, result, std::numeric_limits<number<T, et_on> >::epsilon().backend()); |
405 | else | |
406 | eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value()); | |
92f5a8d4 | 407 | if (eval_get_sign(lim) < 0) |
7c673cae FG |
408 | lim.negate(); |
409 | INSTRUMENT_BACKEND(lim); | |
410 | ||
411 | ui_type k = 1; | |
412 | do | |
413 | { | |
414 | ++k; | |
415 | eval_multiply(pow, t); | |
416 | eval_divide(t2, pow, k); | |
417 | INSTRUMENT_BACKEND(t2); | |
92f5a8d4 | 418 | if (alternate && ((k & 1) != 0)) |
7c673cae FG |
419 | eval_add(result, t2); |
420 | else | |
421 | eval_subtract(result, t2); | |
422 | INSTRUMENT_BACKEND(result); | |
92f5a8d4 | 423 | } while (lim.compare(t2) < 0); |
7c673cae FG |
424 | } |
425 | ||
426 | template <class T> | |
427 | const T& get_constant_log10() | |
428 | { | |
92f5a8d4 TL |
429 | static BOOST_MP_THREAD_LOCAL T result; |
430 | static BOOST_MP_THREAD_LOCAL long digits = 0; | |
431 | #ifndef BOOST_MP_USING_THREAD_LOCAL | |
7c673cae | 432 | static BOOST_MP_THREAD_LOCAL bool b = false; |
92f5a8d4 TL |
433 | constant_initializer<T, &get_constant_log10<T> >::do_nothing(); |
434 | ||
435 | if (!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value())) | |
7c673cae | 436 | { |
92f5a8d4 TL |
437 | b = true; |
438 | #else | |
439 | if ((digits != boost::multiprecision::detail::digits2<number<T> >::value())) | |
440 | { | |
441 | #endif | |
7c673cae | 442 | typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; |
92f5a8d4 | 443 | T ten; |
7c673cae FG |
444 | ten = ui_type(10u); |
445 | eval_log(result, ten); | |
7c673cae FG |
446 | digits = boost::multiprecision::detail::digits2<number<T> >::value(); |
447 | } | |
448 | ||
7c673cae FG |
449 | return result; |
450 | } | |
451 | ||
452 | template <class T> | |
453 | void eval_log10(T& result, const T& arg) | |
454 | { | |
455 | BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log10 function is only valid for floating point types."); | |
456 | eval_log(result, arg); | |
457 | eval_divide(result, get_constant_log10<T>()); | |
458 | } | |
459 | ||
460 | template <class R, class T> | |
461 | inline void eval_log2(R& result, const T& a) | |
462 | { | |
463 | eval_log(result, a); | |
464 | eval_divide(result, get_constant_ln2<R>()); | |
465 | } | |
466 | ||
92f5a8d4 | 467 | template <typename T> |
7c673cae FG |
468 | inline void eval_pow(T& result, const T& x, const T& a) |
469 | { | |
470 | BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The pow function is only valid for floating point types."); | |
471 | typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type; | |
92f5a8d4 | 472 | typedef typename mpl::front<typename T::float_types>::type fp_type; |
7c673cae | 473 | |
92f5a8d4 | 474 | if ((&result == &x) || (&result == &a)) |
7c673cae FG |
475 | { |
476 | T t; | |
477 | eval_pow(t, x, a); | |
478 | result = t; | |
479 | return; | |
480 | } | |
481 | ||
92f5a8d4 | 482 | if ((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0)) |
7c673cae FG |
483 | { |
484 | result = x; | |
485 | return; | |
486 | } | |
92f5a8d4 | 487 | if (a.compare(si_type(0)) == 0) |
b32b8144 FG |
488 | { |
489 | result = si_type(1); | |
490 | return; | |
491 | } | |
7c673cae FG |
492 | |
493 | int type = eval_fpclassify(x); | |
494 | ||
92f5a8d4 | 495 | switch (type) |
7c673cae | 496 | { |
7c673cae | 497 | case FP_ZERO: |
92f5a8d4 | 498 | switch (eval_fpclassify(a)) |
7c673cae FG |
499 | { |
500 | case FP_ZERO: | |
501 | result = si_type(1); | |
502 | break; | |
503 | case FP_NAN: | |
504 | result = a; | |
505 | break; | |
b32b8144 FG |
506 | case FP_NORMAL: |
507 | { | |
508 | // Need to check for a an odd integer as a special case: | |
92f5a8d4 | 509 | try |
b32b8144 FG |
510 | { |
511 | typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type i; | |
512 | eval_convert_to(&i, a); | |
92f5a8d4 | 513 | if (a.compare(i) == 0) |
b32b8144 | 514 | { |
92f5a8d4 | 515 | if (eval_signbit(a)) |
b32b8144 | 516 | { |
92f5a8d4 | 517 | if (i & 1) |
b32b8144 FG |
518 | { |
519 | result = std::numeric_limits<number<T> >::infinity().backend(); | |
92f5a8d4 | 520 | if (eval_signbit(x)) |
b32b8144 FG |
521 | result.negate(); |
522 | errno = ERANGE; | |
523 | } | |
524 | else | |
525 | { | |
526 | result = std::numeric_limits<number<T> >::infinity().backend(); | |
92f5a8d4 | 527 | errno = ERANGE; |
b32b8144 FG |
528 | } |
529 | } | |
92f5a8d4 | 530 | else if (i & 1) |
b32b8144 FG |
531 | { |
532 | result = x; | |
533 | } | |
534 | else | |
535 | result = si_type(0); | |
536 | return; | |
537 | } | |
538 | } | |
92f5a8d4 | 539 | catch (const std::exception&) |
b32b8144 FG |
540 | { |
541 | // fallthrough.. | |
542 | } | |
92f5a8d4 | 543 | BOOST_FALLTHROUGH; |
b32b8144 | 544 | } |
7c673cae | 545 | default: |
92f5a8d4 | 546 | if (eval_signbit(a)) |
b32b8144 FG |
547 | { |
548 | result = std::numeric_limits<number<T> >::infinity().backend(); | |
92f5a8d4 | 549 | errno = ERANGE; |
b32b8144 FG |
550 | } |
551 | else | |
552 | result = x; | |
7c673cae FG |
553 | break; |
554 | } | |
555 | return; | |
556 | case FP_NAN: | |
557 | result = x; | |
92f5a8d4 | 558 | errno = ERANGE; |
7c673cae | 559 | return; |
92f5a8d4 | 560 | default:; |
7c673cae FG |
561 | } |
562 | ||
563 | int s = eval_get_sign(a); | |
92f5a8d4 | 564 | if (s == 0) |
7c673cae FG |
565 | { |
566 | result = si_type(1); | |
567 | return; | |
568 | } | |
569 | ||
92f5a8d4 | 570 | if (s < 0) |
7c673cae FG |
571 | { |
572 | T t, da; | |
573 | t = a; | |
574 | t.negate(); | |
575 | eval_pow(da, x, t); | |
576 | eval_divide(result, si_type(1), da); | |
577 | return; | |
578 | } | |
92f5a8d4 | 579 | |
7c673cae | 580 | typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type an; |
b32b8144 | 581 | typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type max_an = |
92f5a8d4 TL |
582 | std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ? (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::max)() : static_cast<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>(1) << (sizeof(typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type) * CHAR_BIT - 2); |
583 | typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type min_an = | |
584 | std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ? (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::min)() : -min_an; | |
b32b8144 | 585 | |
7c673cae FG |
586 | T fa; |
587 | #ifndef BOOST_NO_EXCEPTIONS | |
588 | try | |
589 | { | |
590 | #endif | |
591 | eval_convert_to(&an, a); | |
92f5a8d4 | 592 | if (a.compare(an) == 0) |
7c673cae FG |
593 | { |
594 | detail::pow_imp(result, x, an, mpl::true_()); | |
595 | return; | |
596 | } | |
597 | #ifndef BOOST_NO_EXCEPTIONS | |
598 | } | |
92f5a8d4 | 599 | catch (const std::exception&) |
7c673cae FG |
600 | { |
601 | // conversion failed, just fall through, value is not an integer. | |
602 | an = (std::numeric_limits<boost::intmax_t>::max)(); | |
603 | } | |
604 | #endif | |
92f5a8d4 | 605 | if ((eval_get_sign(x) < 0)) |
7c673cae FG |
606 | { |
607 | typename boost::multiprecision::detail::canonical<boost::uintmax_t, T>::type aun; | |
608 | #ifndef BOOST_NO_EXCEPTIONS | |
609 | try | |
610 | { | |
611 | #endif | |
612 | eval_convert_to(&aun, a); | |
92f5a8d4 | 613 | if (a.compare(aun) == 0) |
7c673cae FG |
614 | { |
615 | fa = x; | |
616 | fa.negate(); | |
617 | eval_pow(result, fa, a); | |
92f5a8d4 | 618 | if (aun & 1u) |
7c673cae FG |
619 | result.negate(); |
620 | return; | |
621 | } | |
622 | #ifndef BOOST_NO_EXCEPTIONS | |
623 | } | |
92f5a8d4 | 624 | catch (const std::exception&) |
7c673cae FG |
625 | { |
626 | // conversion failed, just fall through, value is not an integer. | |
627 | } | |
628 | #endif | |
b32b8144 FG |
629 | eval_floor(result, a); |
630 | // -1^INF is a special case in C99: | |
92f5a8d4 | 631 | if ((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE)) |
b32b8144 FG |
632 | { |
633 | result = si_type(1); | |
634 | } | |
92f5a8d4 | 635 | else if (a.compare(result) == 0) |
b32b8144 FG |
636 | { |
637 | // exponent is so large we have no fractional part: | |
92f5a8d4 | 638 | if (x.compare(si_type(-1)) < 0) |
b32b8144 FG |
639 | { |
640 | result = std::numeric_limits<number<T, et_on> >::infinity().backend(); | |
641 | } | |
642 | else | |
643 | { | |
644 | result = si_type(0); | |
645 | } | |
646 | } | |
92f5a8d4 | 647 | else if (type == FP_INFINITE) |
b32b8144 FG |
648 | { |
649 | result = std::numeric_limits<number<T, et_on> >::infinity().backend(); | |
650 | } | |
92f5a8d4 | 651 | else if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN) |
b32b8144 | 652 | { |
7c673cae | 653 | result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend(); |
92f5a8d4 | 654 | errno = EDOM; |
b32b8144 | 655 | } |
7c673cae FG |
656 | else |
657 | { | |
658 | BOOST_THROW_EXCEPTION(std::domain_error("Result of pow is undefined or non-real and there is no NaN for this number type.")); | |
659 | } | |
660 | return; | |
661 | } | |
662 | ||
663 | T t, da; | |
664 | ||
665 | eval_subtract(da, a, an); | |
666 | ||
92f5a8d4 | 667 | if ((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an)) |
7c673cae | 668 | { |
92f5a8d4 | 669 | if (a.compare(fp_type(1e-5f)) <= 0) |
7c673cae FG |
670 | { |
671 | // Series expansion for small a. | |
672 | eval_log(t, x); | |
673 | eval_multiply(t, a); | |
674 | hyp0F0(result, t); | |
675 | return; | |
676 | } | |
677 | else | |
678 | { | |
679 | // Series expansion for moderately sized x. Note that for large power of a, | |
680 | // the power of the integer part of a is calculated using the pown function. | |
92f5a8d4 | 681 | if (an) |
7c673cae FG |
682 | { |
683 | da.negate(); | |
684 | t = si_type(1); | |
685 | eval_subtract(t, x); | |
686 | hyp1F0(result, da, t); | |
687 | detail::pow_imp(t, x, an, mpl::true_()); | |
688 | eval_multiply(result, t); | |
689 | } | |
690 | else | |
691 | { | |
692 | da = a; | |
693 | da.negate(); | |
694 | t = si_type(1); | |
695 | eval_subtract(t, x); | |
696 | hyp1F0(result, da, t); | |
697 | } | |
698 | } | |
699 | } | |
700 | else | |
701 | { | |
702 | // Series expansion for pow(x, a). Note that for large power of a, the power | |
703 | // of the integer part of a is calculated using the pown function. | |
92f5a8d4 | 704 | if (an) |
7c673cae FG |
705 | { |
706 | eval_log(t, x); | |
707 | eval_multiply(t, da); | |
708 | eval_exp(result, t); | |
709 | detail::pow_imp(t, x, an, mpl::true_()); | |
710 | eval_multiply(result, t); | |
711 | } | |
712 | else | |
713 | { | |
714 | eval_log(t, x); | |
715 | eval_multiply(t, a); | |
716 | eval_exp(result, t); | |
717 | } | |
718 | } | |
719 | } | |
720 | ||
92f5a8d4 TL |
721 | template <class T, class A> |
722 | #if BOOST_WORKAROUND(BOOST_MSVC, < 1800) | |
723 | inline typename enable_if_c<!is_integral<A>::value, void>::type | |
724 | #else | |
725 | inline typename enable_if_c<is_compatible_arithmetic_type<A, number<T> >::value && !is_integral<A>::value, void>::type | |
726 | #endif | |
727 | eval_pow(T& result, const T& x, const A& a) | |
7c673cae FG |
728 | { |
729 | // Note this one is restricted to float arguments since pow.hpp already has a version for | |
730 | // integer powers.... | |
92f5a8d4 | 731 | typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type; |
7c673cae | 732 | typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type; |
92f5a8d4 | 733 | cast_type c; |
7c673cae FG |
734 | c = a; |
735 | eval_pow(result, x, c); | |
736 | } | |
737 | ||
92f5a8d4 TL |
738 | template <class T, class A> |
739 | #if BOOST_WORKAROUND(BOOST_MSVC, < 1800) | |
740 | inline void | |
741 | #else | |
742 | inline typename enable_if_c<is_compatible_arithmetic_type<A, number<T> >::value, void>::type | |
743 | #endif | |
744 | eval_pow(T& result, const A& x, const T& a) | |
7c673cae | 745 | { |
92f5a8d4 | 746 | typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type; |
7c673cae | 747 | typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type; |
92f5a8d4 | 748 | cast_type c; |
7c673cae FG |
749 | c = x; |
750 | eval_pow(result, c, a); | |
751 | } | |
752 | ||
753 | template <class T> | |
754 | void eval_exp2(T& result, const T& arg) | |
755 | { | |
756 | BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types."); | |
757 | ||
758 | // Check for pure-integer arguments which can be either signed or unsigned. | |
759 | typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i; | |
92f5a8d4 TL |
760 | T temp; |
761 | try | |
762 | { | |
b32b8144 FG |
763 | eval_trunc(temp, arg); |
764 | eval_convert_to(&i, temp); | |
92f5a8d4 | 765 | if (arg.compare(i) == 0) |
b32b8144 FG |
766 | { |
767 | temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u); | |
768 | eval_ldexp(result, temp, i); | |
769 | return; | |
770 | } | |
7c673cae | 771 | } |
92f5a8d4 TL |
772 | catch (const boost::math::rounding_error&) |
773 | { /* Fallthrough */ | |
774 | } | |
775 | catch (const std::runtime_error&) | |
776 | { /* Fallthrough */ | |
777 | } | |
7c673cae FG |
778 | |
779 | temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(2u); | |
780 | eval_pow(result, temp, arg); | |
781 | } | |
782 | ||
92f5a8d4 TL |
783 | namespace detail { |
784 | ||
785 | template <class T> | |
786 | void small_sinh_series(T x, T& result) | |
787 | { | |
788 | typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; | |
789 | bool neg = eval_get_sign(x) < 0; | |
790 | if (neg) | |
791 | x.negate(); | |
792 | T p(x); | |
793 | T mult(x); | |
794 | eval_multiply(mult, x); | |
795 | result = x; | |
796 | ui_type k = 1; | |
797 | ||
798 | T lim(x); | |
799 | eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value()); | |
7c673cae | 800 | |
92f5a8d4 | 801 | do |
7c673cae | 802 | { |
92f5a8d4 TL |
803 | eval_multiply(p, mult); |
804 | eval_divide(p, ++k); | |
805 | eval_divide(p, ++k); | |
806 | eval_add(result, p); | |
807 | } while (p.compare(lim) >= 0); | |
808 | if (neg) | |
809 | result.negate(); | |
810 | } | |
7c673cae | 811 | |
92f5a8d4 TL |
812 | template <class T> |
813 | void sinhcosh(const T& x, T* p_sinh, T* p_cosh) | |
814 | { | |
815 | typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type; | |
816 | typedef typename mpl::front<typename T::float_types>::type fp_type; | |
7c673cae | 817 | |
92f5a8d4 TL |
818 | switch (eval_fpclassify(x)) |
819 | { | |
820 | case FP_NAN: | |
821 | errno = EDOM; | |
822 | // fallthrough... | |
823 | case FP_INFINITE: | |
824 | if (p_sinh) | |
825 | *p_sinh = x; | |
826 | if (p_cosh) | |
7c673cae | 827 | { |
92f5a8d4 TL |
828 | *p_cosh = x; |
829 | if (eval_get_sign(x) < 0) | |
830 | p_cosh->negate(); | |
831 | } | |
832 | return; | |
833 | case FP_ZERO: | |
834 | if (p_sinh) | |
835 | *p_sinh = x; | |
836 | if (p_cosh) | |
837 | *p_cosh = ui_type(1); | |
838 | return; | |
839 | default:; | |
7c673cae FG |
840 | } |
841 | ||
92f5a8d4 TL |
842 | bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0; |
843 | ||
844 | if (p_cosh || !small_sinh) | |
7c673cae | 845 | { |
92f5a8d4 TL |
846 | T e_px, e_mx; |
847 | eval_exp(e_px, x); | |
848 | eval_divide(e_mx, ui_type(1), e_px); | |
849 | if (eval_signbit(e_mx) != eval_signbit(e_px)) | |
850 | e_mx.negate(); // Handles lack of signed zero in some types | |
7c673cae | 851 | |
92f5a8d4 | 852 | if (p_sinh) |
7c673cae | 853 | { |
92f5a8d4 | 854 | if (small_sinh) |
7c673cae | 855 | { |
92f5a8d4 | 856 | small_sinh_series(x, *p_sinh); |
7c673cae | 857 | } |
92f5a8d4 TL |
858 | else |
859 | { | |
860 | eval_subtract(*p_sinh, e_px, e_mx); | |
861 | eval_ldexp(*p_sinh, *p_sinh, -1); | |
7c673cae FG |
862 | } |
863 | } | |
92f5a8d4 | 864 | if (p_cosh) |
7c673cae | 865 | { |
92f5a8d4 TL |
866 | eval_add(*p_cosh, e_px, e_mx); |
867 | eval_ldexp(*p_cosh, *p_cosh, -1); | |
7c673cae FG |
868 | } |
869 | } | |
92f5a8d4 TL |
870 | else |
871 | { | |
872 | small_sinh_series(x, *p_sinh); | |
873 | } | |
874 | } | |
7c673cae FG |
875 | |
876 | } // namespace detail | |
877 | ||
878 | template <class T> | |
879 | inline void eval_sinh(T& result, const T& x) | |
880 | { | |
881 | BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sinh function is only valid for floating point types."); | |
882 | detail::sinhcosh(x, &result, static_cast<T*>(0)); | |
883 | } | |
884 | ||
885 | template <class T> | |
886 | inline void eval_cosh(T& result, const T& x) | |
887 | { | |
888 | BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cosh function is only valid for floating point types."); | |
889 | detail::sinhcosh(x, static_cast<T*>(0), &result); | |
890 | } | |
891 | ||
892 | template <class T> | |
893 | inline void eval_tanh(T& result, const T& x) | |
894 | { | |
895 | BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types."); | |
92f5a8d4 TL |
896 | T c; |
897 | detail::sinhcosh(x, &result, &c); | |
898 | if ((eval_fpclassify(result) == FP_INFINITE) && (eval_fpclassify(c) == FP_INFINITE)) | |
899 | { | |
900 | bool s = eval_signbit(result) != eval_signbit(c); | |
901 | result = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u); | |
902 | if (s) | |
903 | result.negate(); | |
904 | return; | |
905 | } | |
906 | eval_divide(result, c); | |
7c673cae FG |
907 | } |
908 | ||
909 | #ifdef BOOST_MSVC | |
910 | #pragma warning(pop) | |
911 | #endif |