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
13 //
14
15 #ifdef BOOST_MSVC
16 #pragma warning(push)
17 #pragma warning(disable : 6326) // comparison of two constants
18 #endif
19
20 namespace detail {
21
22 template <typename T, typename U>
pow_imp(T & result,const T & t,const U & p,const mpl::false_ &)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
33 if (&result == &t)
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.
42 if (U(p % U(2)) != U(0))
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
54 while (U(p2 /= 2) != U(0))
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
61 if (has_binary_power)
62 {
63 // Multiply the result with each binary power contained in the exponent.
64 eval_multiply(result, x);
65 }
66 }
67 }
68
69 template <typename T, typename U>
pow_imp(T & result,const T & t,const U & p,const mpl::true_ &)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:
73 typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type;
74 typedef typename make_unsigned<U>::type ui_type;
75
76 if (p < 0)
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
90 template <typename T, typename U>
eval_pow(T & result,const T & t,const U & p)91 inline typename enable_if_c<is_integral<U>::value>::type eval_pow(T& result, const T& t, const U& p)
92 {
93 detail::pow_imp(result, t, p, boost::is_signed<U>());
94 }
95
96 template <class T>
hyp0F0(T & H0F0,const T & x)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();
107 T t;
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);
115 if (eval_get_sign(lim) < 0)
116 lim.negate();
117
118 ui_type n;
119
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();
124 // Series expansion of hyperg_0f0(; ; x).
125 for (n = 2; n < series_limit; ++n)
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;
131 if (neg)
132 x_pow_n_div_n_fact.negate();
133 if (lim.compare(x_pow_n_div_n_fact) > 0)
134 break;
135 if (neg)
136 x_pow_n_div_n_fact.negate();
137 }
138 if (n >= series_limit)
139 BOOST_THROW_EXCEPTION(std::runtime_error("H0F0 failed to converge"));
140 }
141
142 template <class T>
hyp1F0(T & H1F0,const T & a,const T & x)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);
156 T pochham_a(a);
157 T ap(a);
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());
163 if (eval_get_sign(lim) < 0)
164 lim.negate();
165
166 si_type n;
167 T term, part;
168
169 const si_type series_limit =
170 boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
171 ? 100
172 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
173 // Series expansion of hyperg_1f0(a; ; x).
174 for (n = 2; n < series_limit; n++)
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);
182 if (eval_get_sign(term) < 0)
183 term.negate();
184 if (lim.compare(term) >= 0)
185 break;
186 }
187 if (n >= series_limit)
188 BOOST_THROW_EXCEPTION(std::runtime_error("H1F0 failed to converge"));
189 }
190
191 template <class T>
eval_exp(T & result,const T & x)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.");
195 if (&x == &result)
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;
203 typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
204 typedef typename T::exponent_type exp_type;
205 typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;
206
207 // Handle special arguments.
208 int type = eval_fpclassify(x);
209 bool isneg = eval_get_sign(x) < 0;
210 if (type == (int)FP_NAN)
211 {
212 result = x;
213 errno = EDOM;
214 return;
215 }
216 else if (type == (int)FP_INFINITE)
217 {
218 if (isneg)
219 result = ui_type(0u);
220 else
221 result = x;
222 return;
223 }
224 else if (type == (int)FP_ZERO)
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;
233 if (isneg)
234 xx.negate();
235
236 // Check the range of the argument.
237 if (xx.compare(si_type(1)) <= 0)
238 {
239 //
240 // Use series for exp(x) - 1:
241 //
242 T lim;
243 if (std::numeric_limits<number<T, et_on> >::is_specialized)
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;
252 result = si_type(1);
253 if (isneg)
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);
260 while (exp_series.compare(lim) > 0)
261 {
262 ++k;
263 eval_multiply(exp_series, xx);
264 eval_divide(exp_series, ui_type(k));
265 if (isneg && (k & 1))
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);
277 if (x.compare(ll) == 0)
278 {
279 detail::pow_imp(result, get_constant_e<T>(), ll, mpl::true_());
280 return;
281 }
282 else if (exp_series.compare(x) == 0)
283 {
284 // We have a value that has no fractional part, but is too large to fit
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:
287 if (isneg)
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 }
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
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
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
332 if (isneg)
333 eval_divide(result, ui_type(1), exp_series);
334 else
335 result = exp_series;
336 }
337
338 template <class T>
eval_log(T & result,const T & arg)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;
349 typedef typename T::exponent_type exp_type;
350 typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;
351 typedef typename mpl::front<typename T::float_types>::type fp_type;
352 int s = eval_signbit(arg);
353 switch (eval_fpclassify(arg))
354 {
355 case FP_NAN:
356 result = arg;
357 errno = EDOM;
358 return;
359 case FP_INFINITE:
360 if (s)
361 break;
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 }
370 if (s)
371 {
372 result = std::numeric_limits<number<T> >::quiet_NaN().backend();
373 errno = EDOM;
374 return;
375 }
376
377 exp_type e;
378 T t;
379 eval_frexp(t, arg, &e);
380 bool alternate = false;
381
382 if (t.compare(fp_type(2) / fp_type(3)) <= 0)
383 {
384 alternate = true;
385 eval_ldexp(t, t, 1);
386 --e;
387 }
388
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 */
392 if (!alternate)
393 t.negate(); /* 0 <= t <= 0.33333 */
394 T pow = t;
395 T lim;
396 T t2;
397
398 if (alternate)
399 eval_add(result, t);
400 else
401 eval_subtract(result, t);
402
403 if (std::numeric_limits<number<T, et_on> >::is_specialized)
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());
407 if (eval_get_sign(lim) < 0)
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);
418 if (alternate && ((k & 1) != 0))
419 eval_add(result, t2);
420 else
421 eval_subtract(result, t2);
422 INSTRUMENT_BACKEND(result);
423 } while (lim.compare(t2) < 0);
424 }
425
426 template <class T>
get_constant_log10()427 const T& get_constant_log10()
428 {
429 static BOOST_MP_THREAD_LOCAL T result;
430 static BOOST_MP_THREAD_LOCAL long digits = 0;
431 #ifndef BOOST_MP_USING_THREAD_LOCAL
432 static BOOST_MP_THREAD_LOCAL bool b = false;
433 constant_initializer<T, &get_constant_log10<T> >::do_nothing();
434
435 if (!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
436 {
437 b = true;
438 #else
439 if ((digits != boost::multiprecision::detail::digits2<number<T> >::value()))
440 {
441 #endif
442 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
443 T ten;
444 ten = ui_type(10u);
445 eval_log(result, ten);
446 digits = boost::multiprecision::detail::digits2<number<T> >::value();
447 }
448
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
467 template <typename T>
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;
472 typedef typename mpl::front<typename T::float_types>::type fp_type;
473
474 if ((&result == &x) || (&result == &a))
475 {
476 T t;
477 eval_pow(t, x, a);
478 result = t;
479 return;
480 }
481
482 if ((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0))
483 {
484 result = x;
485 return;
486 }
487 if (a.compare(si_type(0)) == 0)
488 {
489 result = si_type(1);
490 return;
491 }
492
493 int type = eval_fpclassify(x);
494
495 switch (type)
496 {
497 case FP_ZERO:
498 switch (eval_fpclassify(a))
499 {
500 case FP_ZERO:
501 result = si_type(1);
502 break;
503 case FP_NAN:
504 result = a;
505 break;
506 case FP_NORMAL:
507 {
508 // Need to check for a an odd integer as a special case:
509 try
510 {
511 typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type i;
512 eval_convert_to(&i, a);
513 if (a.compare(i) == 0)
514 {
515 if (eval_signbit(a))
516 {
517 if (i & 1)
518 {
519 result = std::numeric_limits<number<T> >::infinity().backend();
520 if (eval_signbit(x))
521 result.negate();
522 errno = ERANGE;
523 }
524 else
525 {
526 result = std::numeric_limits<number<T> >::infinity().backend();
527 errno = ERANGE;
528 }
529 }
530 else if (i & 1)
531 {
532 result = x;
533 }
534 else
535 result = si_type(0);
536 return;
537 }
538 }
539 catch (const std::exception&)
540 {
541 // fallthrough..
542 }
543 BOOST_FALLTHROUGH;
544 }
545 default:
546 if (eval_signbit(a))
547 {
548 result = std::numeric_limits<number<T> >::infinity().backend();
549 errno = ERANGE;
550 }
551 else
552 result = x;
553 break;
554 }
555 return;
556 case FP_NAN:
557 result = x;
558 errno = ERANGE;
559 return;
560 default:;
561 }
562
563 int s = eval_get_sign(a);
564 if (s == 0)
565 {
566 result = si_type(1);
567 return;
568 }
569
570 if (s < 0)
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 }
579
580 typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type an;
581 typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type max_an =
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;
585
586 T fa;
587 #ifndef BOOST_NO_EXCEPTIONS
588 try
589 {
590 #endif
591 eval_convert_to(&an, a);
592 if (a.compare(an) == 0)
593 {
594 detail::pow_imp(result, x, an, mpl::true_());
595 return;
596 }
597 #ifndef BOOST_NO_EXCEPTIONS
598 }
599 catch (const std::exception&)
600 {
601 // conversion failed, just fall through, value is not an integer.
602 an = (std::numeric_limits<boost::intmax_t>::max)();
603 }
604 #endif
605 if ((eval_get_sign(x) < 0))
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);
613 if (a.compare(aun) == 0)
614 {
615 fa = x;
616 fa.negate();
617 eval_pow(result, fa, a);
618 if (aun & 1u)
619 result.negate();
620 return;
621 }
622 #ifndef BOOST_NO_EXCEPTIONS
623 }
624 catch (const std::exception&)
625 {
626 // conversion failed, just fall through, value is not an integer.
627 }
628 #endif
629 eval_floor(result, a);
630 // -1^INF is a special case in C99:
631 if ((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE))
632 {
633 result = si_type(1);
634 }
635 else if (a.compare(result) == 0)
636 {
637 // exponent is so large we have no fractional part:
638 if (x.compare(si_type(-1)) < 0)
639 {
640 result = std::numeric_limits<number<T, et_on> >::infinity().backend();
641 }
642 else
643 {
644 result = si_type(0);
645 }
646 }
647 else if (type == FP_INFINITE)
648 {
649 result = std::numeric_limits<number<T, et_on> >::infinity().backend();
650 }
651 else if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
652 {
653 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
654 errno = EDOM;
655 }
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
667 if ((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an))
668 {
669 if (a.compare(fp_type(1e-5f)) <= 0)
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.
681 if (an)
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.
704 if (an)
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
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)
728 {
729 // Note this one is restricted to float arguments since pow.hpp already has a version for
730 // integer powers....
731 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
732 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
733 cast_type c;
734 c = a;
735 eval_pow(result, x, c);
736 }
737
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)
745 {
746 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
747 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
748 cast_type c;
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;
760 T temp;
761 try
762 {
763 eval_trunc(temp, arg);
764 eval_convert_to(&i, temp);
765 if (arg.compare(i) == 0)
766 {
767 temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
768 eval_ldexp(result, temp, i);
769 return;
770 }
771 }
772 catch (const boost::math::rounding_error&)
773 { /* Fallthrough */
774 }
775 catch (const std::runtime_error&)
776 { /* Fallthrough */
777 }
778
779 temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(2u);
780 eval_pow(result, temp, arg);
781 }
782
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());
800
801 do
802 {
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 }
811
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;
817
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)
827 {
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:;
840 }
841
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)
845 {
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
851
852 if (p_sinh)
853 {
854 if (small_sinh)
855 {
856 small_sinh_series(x, *p_sinh);
857 }
858 else
859 {
860 eval_subtract(*p_sinh, e_px, e_mx);
861 eval_ldexp(*p_sinh, *p_sinh, -1);
862 }
863 }
864 if (p_cosh)
865 {
866 eval_add(*p_cosh, e_px, e_mx);
867 eval_ldexp(*p_cosh, *p_cosh, -1);
868 }
869 }
870 else
871 {
872 small_sinh_series(x, *p_sinh);
873 }
874 }
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.");
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);
907 }
908
909 #ifdef BOOST_MSVC
910 #pragma warning(pop)
911 #endif
912