]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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> | |
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> | |
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> | |
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> | |
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> | |
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> | |
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; | |
b32b8144 | 211 | errno = EDOM; |
7c673cae FG |
212 | return; |
213 | } | |
214 | else if(type == (int)FP_INFINITE) | |
215 | { | |
7c673cae FG |
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 | } | |
b32b8144 FG |
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 | } | |
7c673cae FG |
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 | ||
b32b8144 FG |
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 | ||
7c673cae FG |
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> | |
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; | |
b32b8144 FG |
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 | } | |
7c673cae FG |
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> | |
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> | |
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> | |
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> | |
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 | ||
b32b8144 | 474 | if((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0)) |
7c673cae FG |
475 | { |
476 | result = x; | |
477 | return; | |
478 | } | |
b32b8144 FG |
479 | if(a.compare(si_type(0)) == 0) |
480 | { | |
481 | result = si_type(1); | |
482 | return; | |
483 | } | |
7c673cae FG |
484 | |
485 | int type = eval_fpclassify(x); | |
486 | ||
487 | switch(type) | |
488 | { | |
7c673cae FG |
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; | |
b32b8144 FG |
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 | } | |
7c673cae | 536 | default: |
b32b8144 FG |
537 | if(eval_signbit(a)) |
538 | { | |
539 | result = std::numeric_limits<number<T> >::infinity().backend(); | |
540 | errno = ERANGE; | |
541 | } | |
542 | else | |
543 | result = x; | |
7c673cae FG |
544 | break; |
545 | } | |
546 | return; | |
547 | case FP_NAN: | |
548 | result = x; | |
b32b8144 | 549 | errno = ERANGE; |
7c673cae FG |
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; | |
b32b8144 FG |
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 | ||
7c673cae FG |
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 | |
b32b8144 FG |
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 | { | |
7c673cae | 649 | result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend(); |
b32b8144 FG |
650 | errno = EDOM; |
651 | } | |
7c673cae FG |
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 | ||
b32b8144 | 663 | if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an)) |
7c673cae FG |
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> | |
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> | |
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> | |
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; | |
b32b8144 FG |
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 | } | |
7c673cae | 756 | } |
b32b8144 FG |
757 | catch(const boost::math::rounding_error&) |
758 | { /* Fallthrough */ } | |
759 | catch(const std::runtime_error&) | |
760 | { /* Fallthrough */ } | |
7c673cae FG |
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> | |
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> | |
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: | |
b32b8144 FG |
804 | errno = EDOM; |
805 | // fallthrough... | |
7c673cae FG |
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); | |
b32b8144 FG |
832 | if(eval_signbit(e_mx) != eval_signbit(e_px)) |
833 | e_mx.negate(); // Handles lack of signed zero in some types | |
7c673cae FG |
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> | |
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> | |
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> | |
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); | |
b32b8144 FG |
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 | } | |
7c673cae FG |
889 | eval_divide(result, c); |
890 | } | |
891 | ||
892 | #ifdef BOOST_MSVC | |
893 | #pragma warning(pop) | |
894 | #endif |