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