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