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

source code of boost/libs/multiprecision/include/boost/multiprecision/detail/functions/pow.hpp