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<is_integral<U> >::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 bool b = false;
428    static BOOST_MP_THREAD_LOCAL long digits = boost::multiprecision::detail::digits2<number<T> >::value();
429    if(!b || (digits != boost::multiprecision::detail::digits2<number<T> >::value()))
430    {
431       typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
432       T ten;
433       ten = ui_type(10u);
434       eval_log(result, ten);
435       b = true;
436       digits = boost::multiprecision::detail::digits2<number<T> >::value();
437    }
438 
439    constant_initializer<T, &get_constant_log10<T> >::do_nothing();
440 
441    return result;
442 }
443 
444 template <class T>
eval_log10(T & result,const T & arg)445 void eval_log10(T& result, const T& arg)
446 {
447    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log10 function is only valid for floating point types.");
448    eval_log(result, arg);
449    eval_divide(result, get_constant_log10<T>());
450 }
451 
452 template <class R, class T>
eval_log2(R & result,const T & a)453 inline void eval_log2(R& result, const T& a)
454 {
455    eval_log(result, a);
456    eval_divide(result, get_constant_ln2<R>());
457 }
458 
459 template<typename T>
eval_pow(T & result,const T & x,const T & a)460 inline void eval_pow(T& result, const T& x, const T& a)
461 {
462    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The pow function is only valid for floating point types.");
463    typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
464    typedef typename mpl::front<typename T::float_types>::type fp_type;
465 
466    if((&result == &x) || (&result == &a))
467    {
468       T t;
469       eval_pow(t, x, a);
470       result = t;
471       return;
472    }
473 
474    if((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0))
475    {
476       result = x;
477       return;
478    }
479    if(a.compare(si_type(0)) == 0)
480    {
481       result = si_type(1);
482       return;
483    }
484 
485    int type = eval_fpclassify(x);
486 
487    switch(type)
488    {
489    case FP_ZERO:
490       switch(eval_fpclassify(a))
491       {
492       case FP_ZERO:
493          result = si_type(1);
494          break;
495       case FP_NAN:
496          result = a;
497          break;
498       case FP_NORMAL:
499       {
500          // Need to check for a an odd integer as a special case:
501          try
502          {
503             typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type i;
504             eval_convert_to(&i, a);
505             if(a.compare(i) == 0)
506             {
507                if(eval_signbit(a))
508                {
509                   if(i & 1)
510                   {
511                      result = std::numeric_limits<number<T> >::infinity().backend();
512                      if(eval_signbit(x))
513                         result.negate();
514                      errno = ERANGE;
515                   }
516                   else
517                   {
518                      result = std::numeric_limits<number<T> >::infinity().backend();
519                      errno = ERANGE;
520                   }
521                }
522                else if(i & 1)
523                {
524                   result = x;
525                }
526                else
527                   result = si_type(0);
528                return;
529             }
530          }
531          catch(const std::exception&)
532          {
533             // fallthrough..
534          }
535       }
536       default:
537          if(eval_signbit(a))
538          {
539             result = std::numeric_limits<number<T> >::infinity().backend();
540             errno = ERANGE;
541          }
542          else
543             result = x;
544          break;
545       }
546       return;
547    case FP_NAN:
548       result = x;
549       errno = ERANGE;
550       return;
551    default: ;
552    }
553 
554    int s = eval_get_sign(a);
555    if(s == 0)
556    {
557       result = si_type(1);
558       return;
559    }
560 
561    if(s < 0)
562    {
563       T t, da;
564       t = a;
565       t.negate();
566       eval_pow(da, x, t);
567       eval_divide(result, si_type(1), da);
568       return;
569    }
570 
571    typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type an;
572    typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type max_an =
573       std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ?
574       (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::max)() :
575       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);
576    typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type min_an =
577       std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ?
578       (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::min)() :
579       -min_an;
580 
581 
582    T fa;
583 #ifndef BOOST_NO_EXCEPTIONS
584    try
585    {
586 #endif
587       eval_convert_to(&an, a);
588       if(a.compare(an) == 0)
589       {
590          detail::pow_imp(result, x, an, mpl::true_());
591          return;
592       }
593 #ifndef BOOST_NO_EXCEPTIONS
594    }
595    catch(const std::exception&)
596    {
597       // conversion failed, just fall through, value is not an integer.
598       an = (std::numeric_limits<boost::intmax_t>::max)();
599    }
600 #endif
601    if((eval_get_sign(x) < 0))
602    {
603       typename boost::multiprecision::detail::canonical<boost::uintmax_t, T>::type aun;
604 #ifndef BOOST_NO_EXCEPTIONS
605       try
606       {
607 #endif
608          eval_convert_to(&aun, a);
609          if(a.compare(aun) == 0)
610          {
611             fa = x;
612             fa.negate();
613             eval_pow(result, fa, a);
614             if(aun & 1u)
615                result.negate();
616             return;
617          }
618 #ifndef BOOST_NO_EXCEPTIONS
619       }
620       catch(const std::exception&)
621       {
622          // conversion failed, just fall through, value is not an integer.
623       }
624 #endif
625       eval_floor(result, a);
626       // -1^INF is a special case in C99:
627       if((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE))
628       {
629          result = si_type(1);
630       }
631       else if(a.compare(result) == 0)
632       {
633          // exponent is so large we have no fractional part:
634          if(x.compare(si_type(-1)) < 0)
635          {
636             result = std::numeric_limits<number<T, et_on> >::infinity().backend();
637          }
638          else
639          {
640             result = si_type(0);
641          }
642       }
643       else if(type == FP_INFINITE)
644       {
645          result = std::numeric_limits<number<T, et_on> >::infinity().backend();
646       }
647       else if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
648       {
649          result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
650          errno = EDOM;
651       }
652       else
653       {
654          BOOST_THROW_EXCEPTION(std::domain_error("Result of pow is undefined or non-real and there is no NaN for this number type."));
655       }
656       return;
657    }
658 
659    T t, da;
660 
661    eval_subtract(da, a, an);
662 
663    if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an))
664    {
665       if(a.compare(fp_type(1e-5f)) <= 0)
666       {
667          // Series expansion for small a.
668          eval_log(t, x);
669          eval_multiply(t, a);
670          hyp0F0(result, t);
671          return;
672       }
673       else
674       {
675          // Series expansion for moderately sized x. Note that for large power of a,
676          // the power of the integer part of a is calculated using the pown function.
677          if(an)
678          {
679             da.negate();
680             t = si_type(1);
681             eval_subtract(t, x);
682             hyp1F0(result, da, t);
683             detail::pow_imp(t, x, an, mpl::true_());
684             eval_multiply(result, t);
685          }
686          else
687          {
688             da = a;
689             da.negate();
690             t = si_type(1);
691             eval_subtract(t, x);
692             hyp1F0(result, da, t);
693          }
694       }
695    }
696    else
697    {
698       // Series expansion for pow(x, a). Note that for large power of a, the power
699       // of the integer part of a is calculated using the pown function.
700       if(an)
701       {
702          eval_log(t, x);
703          eval_multiply(t, da);
704          eval_exp(result, t);
705          detail::pow_imp(t, x, an, mpl::true_());
706          eval_multiply(result, t);
707       }
708       else
709       {
710          eval_log(t, x);
711          eval_multiply(t, a);
712          eval_exp(result, t);
713       }
714    }
715 }
716 
717 template<class T, class A>
eval_pow(T & result,const T & x,const A & a)718 inline typename enable_if<is_floating_point<A>, void>::type eval_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    typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
723    typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
724    cast_type c;
725    c = a;
726    eval_pow(result, x, c);
727 }
728 
729 template<class T, class A>
eval_pow(T & result,const A & x,const T & a)730 inline typename enable_if<is_arithmetic<A>, void>::type eval_pow(T& result, const A& x, const T& a)
731 {
732    typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
733    typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
734    cast_type c;
735    c = x;
736    eval_pow(result, c, a);
737 }
738 
739 template <class T>
eval_exp2(T & result,const T & arg)740 void eval_exp2(T& result, const T& arg)
741 {
742    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The log function is only valid for floating point types.");
743 
744    // Check for pure-integer arguments which can be either signed or unsigned.
745    typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i;
746    T temp;
747    try {
748       eval_trunc(temp, arg);
749       eval_convert_to(&i, temp);
750       if(arg.compare(i) == 0)
751       {
752          temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
753          eval_ldexp(result, temp, i);
754          return;
755       }
756    }
757    catch(const boost::math::rounding_error&)
758    { /* Fallthrough */ }
759    catch(const std::runtime_error&)
760    { /* Fallthrough */ }
761 
762    temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(2u);
763    eval_pow(result, temp, arg);
764 }
765 
766 namespace detail{
767 
768    template <class T>
small_sinh_series(T x,T & result)769    void small_sinh_series(T x, T& result)
770    {
771       typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
772       bool neg = eval_get_sign(x) < 0;
773       if(neg)
774          x.negate();
775       T p(x);
776       T mult(x);
777       eval_multiply(mult, x);
778       result = x;
779       ui_type k = 1;
780 
781       T lim(x);
782       eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
783 
784       do
785       {
786          eval_multiply(p, mult);
787          eval_divide(p, ++k);
788          eval_divide(p, ++k);
789          eval_add(result, p);
790       }while(p.compare(lim) >= 0);
791       if(neg)
792          result.negate();
793    }
794 
795    template <class T>
sinhcosh(const T & x,T * p_sinh,T * p_cosh)796    void sinhcosh(const T& x, T* p_sinh, T* p_cosh)
797    {
798       typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
799       typedef typename mpl::front<typename T::float_types>::type fp_type;
800 
801       switch(eval_fpclassify(x))
802       {
803       case FP_NAN:
804          errno = EDOM;
805          // fallthrough...
806       case FP_INFINITE:
807          if(p_sinh)
808             *p_sinh = x;
809          if(p_cosh)
810          {
811             *p_cosh = x;
812             if(eval_get_sign(x) < 0)
813                p_cosh->negate();
814          }
815          return;
816       case FP_ZERO:
817          if(p_sinh)
818             *p_sinh = x;
819          if(p_cosh)
820             *p_cosh = ui_type(1);
821          return;
822       default: ;
823       }
824 
825       bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0;
826 
827       if(p_cosh || !small_sinh)
828       {
829          T e_px, e_mx;
830          eval_exp(e_px, x);
831          eval_divide(e_mx, ui_type(1), e_px);
832          if(eval_signbit(e_mx) != eval_signbit(e_px))
833             e_mx.negate();  // Handles lack of signed zero in some types
834 
835          if(p_sinh)
836          {
837             if(small_sinh)
838             {
839                small_sinh_series(x, *p_sinh);
840             }
841             else
842             {
843                eval_subtract(*p_sinh, e_px, e_mx);
844                eval_ldexp(*p_sinh, *p_sinh, -1);
845             }
846          }
847          if(p_cosh)
848          {
849             eval_add(*p_cosh, e_px, e_mx);
850             eval_ldexp(*p_cosh, *p_cosh, -1);
851          }
852       }
853       else
854       {
855          small_sinh_series(x, *p_sinh);
856       }
857    }
858 
859 } // namespace detail
860 
861 template <class T>
eval_sinh(T & result,const T & x)862 inline void eval_sinh(T& result, const T& x)
863 {
864    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sinh function is only valid for floating point types.");
865    detail::sinhcosh(x, &result, static_cast<T*>(0));
866 }
867 
868 template <class T>
eval_cosh(T & result,const T & x)869 inline void eval_cosh(T& result, const T& x)
870 {
871    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cosh function is only valid for floating point types.");
872    detail::sinhcosh(x, static_cast<T*>(0), &result);
873 }
874 
875 template <class T>
eval_tanh(T & result,const T & x)876 inline void eval_tanh(T& result, const T& x)
877 {
878    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types.");
879   T c;
880   detail::sinhcosh(x, &result, &c);
881   if((eval_fpclassify(result) == FP_INFINITE) && (eval_fpclassify(c) == FP_INFINITE))
882   {
883      bool s = eval_signbit(result) != eval_signbit(c);
884      result = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
885      if(s)
886         result.negate();
887      return;
888   }
889   eval_divide(result, c);
890 }
891 
892 #ifdef BOOST_MSVC
893 #pragma warning(pop)
894 #endif
895