1
2 // Copyright Christopher Kormanyos 2002 - 2011.
3 // Copyright 2011 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 template <class T>
hyp0F1(T & result,const T & b,const T & x)21 void hyp0F1(T& result, const T& b, const T& x)
22 {
23 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
24 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
25
26 // Compute the series representation of Hypergeometric0F1 taken from
27 // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F1/06/01/01/
28 // There are no checks on input range or parameter boundaries.
29
30 T x_pow_n_div_n_fact(x);
31 T pochham_b(b);
32 T bp(b);
33
34 eval_divide(result, x_pow_n_div_n_fact, pochham_b);
35 eval_add(result, ui_type(1));
36
37 si_type n;
38
39 T tol;
40 tol = ui_type(1);
41 eval_ldexp(tol, tol, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
42 eval_multiply(tol, result);
43 if (eval_get_sign(tol) < 0)
44 tol.negate();
45 T term;
46
47 const int series_limit =
48 boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
49 ? 100
50 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
51 // Series expansion of hyperg_0f1(; b; x).
52 for (n = 2; n < series_limit; ++n)
53 {
54 eval_multiply(x_pow_n_div_n_fact, x);
55 eval_divide(x_pow_n_div_n_fact, n);
56 eval_increment(bp);
57 eval_multiply(pochham_b, bp);
58
59 eval_divide(term, x_pow_n_div_n_fact, pochham_b);
60 eval_add(result, term);
61
62 bool neg_term = eval_get_sign(term) < 0;
63 if (neg_term)
64 term.negate();
65 if (term.compare(tol) <= 0)
66 break;
67 }
68
69 if (n >= series_limit)
70 BOOST_THROW_EXCEPTION(std::runtime_error("H0F1 Failed to Converge"));
71 }
72
73 template <class T, unsigned N, bool b = boost::multiprecision::detail::is_variable_precision<boost::multiprecision::number<T> >::value>
74 struct scoped_N_precision
75 {
76 template <class U>
scoped_N_precisionscoped_N_precision77 scoped_N_precision(U const&) {}
78 template <class U>
reducescoped_N_precision79 void reduce(U&) {}
80 };
81
82 template <class T, unsigned N>
83 struct scoped_N_precision<T, N, true>
84 {
85 unsigned old_precision, old_arg_precision;
scoped_N_precisionscoped_N_precision86 scoped_N_precision(T& arg)
87 {
88 old_precision = T::default_precision();
89 old_arg_precision = arg.precision();
90 T::default_precision(old_arg_precision * N);
91 arg.precision(old_arg_precision * N);
92 }
~scoped_N_precisionscoped_N_precision93 ~scoped_N_precision()
94 {
95 T::default_precision(old_precision);
96 }
reducescoped_N_precision97 void reduce(T& arg)
98 {
99 arg.precision(old_arg_precision);
100 }
101 };
102
103 template <class T>
reduce_n_half_pi(T & arg,const T & n,bool go_down)104 void reduce_n_half_pi(T& arg, const T& n, bool go_down)
105 {
106 //
107 // We need to perform argument reduction at 3 times the precision of arg
108 // in order to ensure a correct result up to arg = 1/epsilon. Beyond that
109 // the value of n will have been incorrectly calculated anyway since it will
110 // have a value greater than 1/epsilon and no longer be an exact integer value.
111 //
112 // More information in ARGUMENT REDUCTION FOR HUGE ARGUMENTS. K C Ng.
113 //
114 // There are two mutually exclusive ways to achieve this, both of which are
115 // supported here:
116 // 1) To define a fixed precision type with 3 times the precision for the calculation.
117 // 2) To dynamically increase the precision of the variables.
118 //
119 typedef typename boost::multiprecision::detail::transcendental_reduction_type<T>::type reduction_type;
120 //
121 // Make a copy of the arg at higher precision:
122 //
123 reduction_type big_arg(arg);
124 //
125 // Dynamically increase precision when supported, this increases the default
126 // and ups the precision of big_arg to match:
127 //
128 scoped_N_precision<T, 3> scoped_precision(big_arg);
129 //
130 // High precision PI:
131 //
132 reduction_type reduction = get_constant_pi<reduction_type>();
133 eval_ldexp(reduction, reduction, -1); // divide by 2
134 eval_multiply(reduction, n);
135 BOOST_MATH_INSTRUMENT_CODE(big_arg.str(10, std::ios_base::scientific));
136 BOOST_MATH_INSTRUMENT_CODE(reduction.str(10, std::ios_base::scientific));
137
138 if (go_down)
139 eval_subtract(big_arg, reduction, big_arg);
140 else
141 eval_subtract(big_arg, reduction);
142 arg = T(big_arg);
143 //
144 // If arg is a variable precision type, then we have just copied the
145 // precision of big_arg s well it's value. Reduce the precision now:
146 //
147 scoped_precision.reduce(arg);
148 BOOST_MATH_INSTRUMENT_CODE(big_arg.str(10, std::ios_base::scientific));
149 BOOST_MATH_INSTRUMENT_CODE(arg.str(10, std::ios_base::scientific));
150 }
151
152 template <class T>
eval_sin(T & result,const T & x)153 void eval_sin(T& result, const T& x)
154 {
155 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sin function is only valid for floating point types.");
156 BOOST_MATH_INSTRUMENT_CODE(x.str(0, std::ios_base::scientific));
157 if (&result == &x)
158 {
159 T temp;
160 eval_sin(temp, x);
161 result = temp;
162 return;
163 }
164
165 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
166 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
167 typedef typename mpl::front<typename T::float_types>::type fp_type;
168
169 switch (eval_fpclassify(x))
170 {
171 case FP_INFINITE:
172 case FP_NAN:
173 if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
174 {
175 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
176 errno = EDOM;
177 }
178 else
179 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
180 return;
181 case FP_ZERO:
182 result = x;
183 return;
184 default:;
185 }
186
187 // Local copy of the argument
188 T xx = x;
189
190 // Analyze and prepare the phase of the argument.
191 // Make a local, positive copy of the argument, xx.
192 // The argument xx will be reduced to 0 <= xx <= pi/2.
193 bool b_negate_sin = false;
194
195 if (eval_get_sign(x) < 0)
196 {
197 xx.negate();
198 b_negate_sin = !b_negate_sin;
199 }
200
201 T n_pi, t;
202 T half_pi = get_constant_pi<T>();
203 eval_ldexp(half_pi, half_pi, -1); // divide by 2
204 // Remove multiples of pi/2.
205 if (xx.compare(half_pi) > 0)
206 {
207 eval_divide(n_pi, xx, half_pi);
208 eval_trunc(n_pi, n_pi);
209 t = ui_type(4);
210 eval_fmod(t, n_pi, t);
211 bool b_go_down = false;
212 if (t.compare(ui_type(1)) == 0)
213 {
214 b_go_down = true;
215 }
216 else if (t.compare(ui_type(2)) == 0)
217 {
218 b_negate_sin = !b_negate_sin;
219 }
220 else if (t.compare(ui_type(3)) == 0)
221 {
222 b_negate_sin = !b_negate_sin;
223 b_go_down = true;
224 }
225
226 if (b_go_down)
227 eval_increment(n_pi);
228 //
229 // If n_pi is > 1/epsilon, then it is no longer an exact integer value
230 // but an approximation. As a result we can no longer reliably reduce
231 // xx to 0 <= xx < pi/2, nor can we tell the sign of the result as we need
232 // n_pi % 4 for that, but that will always be zero in this situation.
233 // We could use a higher precision type for n_pi, along with division at
234 // higher precision, but that's rather expensive. So for now we do not support
235 // this, and will see if anyone complains and has a legitimate use case.
236 //
237 if (n_pi.compare(get_constant_one_over_epsilon<T>()) > 0)
238 {
239 result = ui_type(0);
240 return;
241 }
242
243 reduce_n_half_pi(xx, n_pi, b_go_down);
244 //
245 // Post reduction we may be a few ulp below zero or above pi/2
246 // given that n_pi was calculated at working precision and not
247 // at the higher precision used for reduction. Correct that now:
248 //
249 if (eval_get_sign(xx) < 0)
250 {
251 xx.negate();
252 b_negate_sin = !b_negate_sin;
253 }
254 if (xx.compare(half_pi) > 0)
255 {
256 eval_ldexp(half_pi, half_pi, 1);
257 eval_subtract(xx, half_pi, xx);
258 eval_ldexp(half_pi, half_pi, -1);
259 b_go_down = !b_go_down;
260 }
261
262 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
263 BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
264 BOOST_ASSERT(xx.compare(half_pi) <= 0);
265 BOOST_ASSERT(xx.compare(ui_type(0)) >= 0);
266 }
267
268 t = half_pi;
269 eval_subtract(t, xx);
270
271 const bool b_zero = eval_get_sign(xx) == 0;
272 const bool b_pi_half = eval_get_sign(t) == 0;
273
274 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
275 BOOST_MATH_INSTRUMENT_CODE(t.str(0, std::ios_base::scientific));
276
277 // Check if the reduced argument is very close to 0 or pi/2.
278 const bool b_near_zero = xx.compare(fp_type(1e-1)) < 0;
279 const bool b_near_pi_half = t.compare(fp_type(1e-1)) < 0;
280
281 if (b_zero)
282 {
283 result = ui_type(0);
284 }
285 else if (b_pi_half)
286 {
287 result = ui_type(1);
288 }
289 else if (b_near_zero)
290 {
291 eval_multiply(t, xx, xx);
292 eval_divide(t, si_type(-4));
293 T t2;
294 t2 = fp_type(1.5);
295 hyp0F1(result, t2, t);
296 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
297 eval_multiply(result, xx);
298 }
299 else if (b_near_pi_half)
300 {
301 eval_multiply(t, t);
302 eval_divide(t, si_type(-4));
303 T t2;
304 t2 = fp_type(0.5);
305 hyp0F1(result, t2, t);
306 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
307 }
308 else
309 {
310 // Scale to a small argument for an efficient Taylor series,
311 // implemented as a hypergeometric function. Use a standard
312 // divide by three identity a certain number of times.
313 // Here we use division by 3^9 --> (19683 = 3^9).
314
315 static const si_type n_scale = 9;
316 static const si_type n_three_pow_scale = static_cast<si_type>(19683L);
317
318 eval_divide(xx, n_three_pow_scale);
319
320 // Now with small arguments, we are ready for a series expansion.
321 eval_multiply(t, xx, xx);
322 eval_divide(t, si_type(-4));
323 T t2;
324 t2 = fp_type(1.5);
325 hyp0F1(result, t2, t);
326 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
327 eval_multiply(result, xx);
328
329 // Convert back using multiple angle identity.
330 for (boost::int32_t k = static_cast<boost::int32_t>(0); k < n_scale; k++)
331 {
332 // Rescale the cosine value using the multiple angle identity.
333 eval_multiply(t2, result, ui_type(3));
334 eval_multiply(t, result, result);
335 eval_multiply(t, result);
336 eval_multiply(t, ui_type(4));
337 eval_subtract(result, t2, t);
338 }
339 }
340
341 if (b_negate_sin)
342 result.negate();
343 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
344 }
345
346 template <class T>
eval_cos(T & result,const T & x)347 void eval_cos(T& result, const T& x)
348 {
349 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cos function is only valid for floating point types.");
350 if (&result == &x)
351 {
352 T temp;
353 eval_cos(temp, x);
354 result = temp;
355 return;
356 }
357
358 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
359 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
360
361 switch (eval_fpclassify(x))
362 {
363 case FP_INFINITE:
364 case FP_NAN:
365 if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
366 {
367 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
368 errno = EDOM;
369 }
370 else
371 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
372 return;
373 case FP_ZERO:
374 result = ui_type(1);
375 return;
376 default:;
377 }
378
379 // Local copy of the argument
380 T xx = x;
381
382 // Analyze and prepare the phase of the argument.
383 // Make a local, positive copy of the argument, xx.
384 // The argument xx will be reduced to 0 <= xx <= pi/2.
385 bool b_negate_cos = false;
386
387 if (eval_get_sign(x) < 0)
388 {
389 xx.negate();
390 }
391 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
392
393 T n_pi, t;
394 T half_pi = get_constant_pi<T>();
395 eval_ldexp(half_pi, half_pi, -1); // divide by 2
396 // Remove even multiples of pi.
397 if (xx.compare(half_pi) > 0)
398 {
399 eval_divide(t, xx, half_pi);
400 eval_trunc(n_pi, t);
401 BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
402 t = ui_type(4);
403 eval_fmod(t, n_pi, t);
404
405 bool b_go_down = false;
406 if (t.compare(ui_type(0)) == 0)
407 {
408 b_go_down = true;
409 }
410 else if (t.compare(ui_type(1)) == 0)
411 {
412 b_negate_cos = true;
413 }
414 else if (t.compare(ui_type(2)) == 0)
415 {
416 b_go_down = true;
417 b_negate_cos = true;
418 }
419 else
420 {
421 BOOST_ASSERT(t.compare(ui_type(3)) == 0);
422 }
423
424 if (b_go_down)
425 eval_increment(n_pi);
426 //
427 // If n_pi is > 1/epsilon, then it is no longer an exact integer value
428 // but an approximation. As a result we can no longer reliably reduce
429 // xx to 0 <= xx < pi/2, nor can we tell the sign of the result as we need
430 // n_pi % 4 for that, but that will always be zero in this situation.
431 // We could use a higher precision type for n_pi, along with division at
432 // higher precision, but that's rather expensive. So for now we do not support
433 // this, and will see if anyone complains and has a legitimate use case.
434 //
435 if (n_pi.compare(get_constant_one_over_epsilon<T>()) > 0)
436 {
437 result = ui_type(1);
438 return;
439 }
440
441 reduce_n_half_pi(xx, n_pi, b_go_down);
442 //
443 // Post reduction we may be a few ulp below zero or above pi/2
444 // given that n_pi was calculated at working precision and not
445 // at the higher precision used for reduction. Correct that now:
446 //
447 if (eval_get_sign(xx) < 0)
448 {
449 xx.negate();
450 b_negate_cos = !b_negate_cos;
451 }
452 if (xx.compare(half_pi) > 0)
453 {
454 eval_ldexp(half_pi, half_pi, 1);
455 eval_subtract(xx, half_pi, xx);
456 eval_ldexp(half_pi, half_pi, -1);
457 }
458 BOOST_ASSERT(xx.compare(half_pi) <= 0);
459 BOOST_ASSERT(xx.compare(ui_type(0)) >= 0);
460 }
461 else
462 {
463 n_pi = ui_type(1);
464 reduce_n_half_pi(xx, n_pi, true);
465 }
466
467 const bool b_zero = eval_get_sign(xx) == 0;
468
469 if (b_zero)
470 {
471 result = si_type(0);
472 }
473 else
474 {
475 eval_sin(result, xx);
476 }
477 if (b_negate_cos)
478 result.negate();
479 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
480 }
481
482 template <class T>
eval_tan(T & result,const T & x)483 void eval_tan(T& result, const T& x)
484 {
485 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tan function is only valid for floating point types.");
486 if (&result == &x)
487 {
488 T temp;
489 eval_tan(temp, x);
490 result = temp;
491 return;
492 }
493 T t;
494 eval_sin(result, x);
495 eval_cos(t, x);
496 eval_divide(result, t);
497 }
498
499 template <class T>
hyp2F1(T & result,const T & a,const T & b,const T & c,const T & x)500 void hyp2F1(T& result, const T& a, const T& b, const T& c, const T& x)
501 {
502 // Compute the series representation of hyperg_2f1 taken from
503 // Abramowitz and Stegun 15.1.1.
504 // There are no checks on input range or parameter boundaries.
505
506 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
507
508 T x_pow_n_div_n_fact(x);
509 T pochham_a(a);
510 T pochham_b(b);
511 T pochham_c(c);
512 T ap(a);
513 T bp(b);
514 T cp(c);
515
516 eval_multiply(result, pochham_a, pochham_b);
517 eval_divide(result, pochham_c);
518 eval_multiply(result, x_pow_n_div_n_fact);
519 eval_add(result, ui_type(1));
520
521 T lim;
522 eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
523
524 if (eval_get_sign(lim) < 0)
525 lim.negate();
526
527 ui_type n;
528 T term;
529
530 const unsigned series_limit =
531 boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
532 ? 100
533 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
534 // Series expansion of hyperg_2f1(a, b; c; x).
535 for (n = 2; n < series_limit; ++n)
536 {
537 eval_multiply(x_pow_n_div_n_fact, x);
538 eval_divide(x_pow_n_div_n_fact, n);
539
540 eval_increment(ap);
541 eval_multiply(pochham_a, ap);
542 eval_increment(bp);
543 eval_multiply(pochham_b, bp);
544 eval_increment(cp);
545 eval_multiply(pochham_c, cp);
546
547 eval_multiply(term, pochham_a, pochham_b);
548 eval_divide(term, pochham_c);
549 eval_multiply(term, x_pow_n_div_n_fact);
550 eval_add(result, term);
551
552 if (eval_get_sign(term) < 0)
553 term.negate();
554 if (lim.compare(term) >= 0)
555 break;
556 }
557 if (n > series_limit)
558 BOOST_THROW_EXCEPTION(std::runtime_error("H2F1 failed to converge."));
559 }
560
561 template <class T>
eval_asin(T & result,const T & x)562 void eval_asin(T& result, const T& x)
563 {
564 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The asin function is only valid for floating point types.");
565 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
566 typedef typename mpl::front<typename T::float_types>::type fp_type;
567
568 if (&result == &x)
569 {
570 T t(x);
571 eval_asin(result, t);
572 return;
573 }
574
575 switch (eval_fpclassify(x))
576 {
577 case FP_NAN:
578 case FP_INFINITE:
579 if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
580 {
581 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
582 errno = EDOM;
583 }
584 else
585 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
586 return;
587 case FP_ZERO:
588 result = x;
589 return;
590 default:;
591 }
592
593 const bool b_neg = eval_get_sign(x) < 0;
594
595 T xx(x);
596 if (b_neg)
597 xx.negate();
598
599 int c = xx.compare(ui_type(1));
600 if (c > 0)
601 {
602 if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
603 {
604 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
605 errno = EDOM;
606 }
607 else
608 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
609 return;
610 }
611 else if (c == 0)
612 {
613 result = get_constant_pi<T>();
614 eval_ldexp(result, result, -1);
615 if (b_neg)
616 result.negate();
617 return;
618 }
619
620 if (xx.compare(fp_type(1e-3)) < 0)
621 {
622 // http://functions.wolfram.com/ElementaryFunctions/ArcSin/26/01/01/
623 eval_multiply(xx, xx);
624 T t1, t2;
625 t1 = fp_type(0.5f);
626 t2 = fp_type(1.5f);
627 hyp2F1(result, t1, t1, t2, xx);
628 eval_multiply(result, x);
629 return;
630 }
631 else if (xx.compare(fp_type(1 - 5e-2f)) > 0)
632 {
633 // http://functions.wolfram.com/ElementaryFunctions/ArcSin/26/01/01/
634 // This branch is simlilar in complexity to Newton iterations down to
635 // the above limit. It is *much* more accurate.
636 T dx1;
637 T t1, t2;
638 eval_subtract(dx1, ui_type(1), xx);
639 t1 = fp_type(0.5f);
640 t2 = fp_type(1.5f);
641 eval_ldexp(dx1, dx1, -1);
642 hyp2F1(result, t1, t1, t2, dx1);
643 eval_ldexp(dx1, dx1, 2);
644 eval_sqrt(t1, dx1);
645 eval_multiply(result, t1);
646 eval_ldexp(t1, get_constant_pi<T>(), -1);
647 result.negate();
648 eval_add(result, t1);
649 if (b_neg)
650 result.negate();
651 return;
652 }
653 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
654 typedef typename boost::multiprecision::detail::canonical<long double, T>::type guess_type;
655 #else
656 typedef fp_type guess_type;
657 #endif
658 // Get initial estimate using standard math function asin.
659 guess_type dd;
660 eval_convert_to(&dd, xx);
661
662 result = (guess_type)(std::asin(dd));
663
664 // Newton-Raphson iteration, we should double our precision with each iteration,
665 // in practice this seems to not quite work in all cases... so terminate when we
666 // have at least 2/3 of the digits correct on the assumption that the correction
667 // we've just added will finish the job...
668
669 boost::intmax_t current_precision = eval_ilogb(result);
670 boost::intmax_t target_precision = std::numeric_limits<number<T> >::is_specialized ?
671 current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3
672 : current_precision - 1 - (boost::multiprecision::detail::digits2<number<T> >::value() * 2) / 3;
673
674 // Newton-Raphson iteration
675 while (current_precision > target_precision)
676 {
677 T sine, cosine;
678 eval_sin(sine, result);
679 eval_cos(cosine, result);
680 eval_subtract(sine, xx);
681 eval_divide(sine, cosine);
682 eval_subtract(result, sine);
683 current_precision = eval_ilogb(sine);
684 if (current_precision <= (std::numeric_limits<typename T::exponent_type>::min)() + 1)
685 break;
686 }
687 if (b_neg)
688 result.negate();
689 }
690
691 template <class T>
eval_acos(T & result,const T & x)692 inline void eval_acos(T& result, const T& x)
693 {
694 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The acos function is only valid for floating point types.");
695 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
696
697 switch (eval_fpclassify(x))
698 {
699 case FP_NAN:
700 case FP_INFINITE:
701 if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
702 {
703 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
704 errno = EDOM;
705 }
706 else
707 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
708 return;
709 case FP_ZERO:
710 result = get_constant_pi<T>();
711 eval_ldexp(result, result, -1); // divide by two.
712 return;
713 }
714
715 T xx;
716 eval_abs(xx, x);
717 int c = xx.compare(ui_type(1));
718
719 if (c > 0)
720 {
721 if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
722 {
723 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
724 errno = EDOM;
725 }
726 else
727 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
728 return;
729 }
730 else if (c == 0)
731 {
732 if (eval_get_sign(x) < 0)
733 result = get_constant_pi<T>();
734 else
735 result = ui_type(0);
736 return;
737 }
738
739 typedef typename mpl::front<typename T::float_types>::type fp_type;
740
741 if (xx.compare(fp_type(1e-3)) < 0)
742 {
743 // https://functions.wolfram.com/ElementaryFunctions/ArcCos/26/01/01/
744 eval_multiply(xx, xx);
745 T t1, t2;
746 t1 = fp_type(0.5f);
747 t2 = fp_type(1.5f);
748 hyp2F1(result, t1, t1, t2, xx);
749 eval_multiply(result, x);
750 eval_ldexp(t1, get_constant_pi<T>(), -1);
751 result.negate();
752 eval_add(result, t1);
753 return;
754 }
755 if (eval_get_sign(x) < 0)
756 {
757 eval_acos(result, xx);
758 result.negate();
759 eval_add(result, get_constant_pi<T>());
760 return;
761 }
762 else if (xx.compare(fp_type(0.85)) > 0)
763 {
764 // https://functions.wolfram.com/ElementaryFunctions/ArcCos/26/01/01/
765 // This branch is simlilar in complexity to Newton iterations down to
766 // the above limit. It is *much* more accurate.
767 T dx1;
768 T t1, t2;
769 eval_subtract(dx1, ui_type(1), xx);
770 t1 = fp_type(0.5f);
771 t2 = fp_type(1.5f);
772 eval_ldexp(dx1, dx1, -1);
773 hyp2F1(result, t1, t1, t2, dx1);
774 eval_ldexp(dx1, dx1, 2);
775 eval_sqrt(t1, dx1);
776 eval_multiply(result, t1);
777 return;
778 }
779
780 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
781 typedef typename boost::multiprecision::detail::canonical<long double, T>::type guess_type;
782 #else
783 typedef fp_type guess_type;
784 #endif
785 // Get initial estimate using standard math function asin.
786 guess_type dd;
787 eval_convert_to(&dd, xx);
788
789 result = (guess_type)(std::acos(dd));
790
791 // Newton-Raphson iteration, we should double our precision with each iteration,
792 // in practice this seems to not quite work in all cases... so terminate when we
793 // have at least 2/3 of the digits correct on the assumption that the correction
794 // we've just added will finish the job...
795
796 boost::intmax_t current_precision = eval_ilogb(result);
797 boost::intmax_t target_precision = std::numeric_limits<number<T> >::is_specialized ?
798 current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3
799 : current_precision - 1 - (boost::multiprecision::detail::digits2<number<T> >::value() * 2) / 3;
800
801 // Newton-Raphson iteration
802 while (current_precision > target_precision)
803 {
804 T sine, cosine;
805 eval_sin(sine, result);
806 eval_cos(cosine, result);
807 eval_subtract(cosine, xx);
808 cosine.negate();
809 eval_divide(cosine, sine);
810 eval_subtract(result, cosine);
811 current_precision = eval_ilogb(cosine);
812 if (current_precision <= (std::numeric_limits<typename T::exponent_type>::min)() + 1)
813 break;
814 }
815 }
816
817 template <class T>
eval_atan(T & result,const T & x)818 void eval_atan(T& result, const T& x)
819 {
820 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan function is only valid for floating point types.");
821 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
822 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
823 typedef typename mpl::front<typename T::float_types>::type fp_type;
824
825 switch (eval_fpclassify(x))
826 {
827 case FP_NAN:
828 result = x;
829 errno = EDOM;
830 return;
831 case FP_ZERO:
832 result = x;
833 return;
834 case FP_INFINITE:
835 if (eval_get_sign(x) < 0)
836 {
837 eval_ldexp(result, get_constant_pi<T>(), -1);
838 result.negate();
839 }
840 else
841 eval_ldexp(result, get_constant_pi<T>(), -1);
842 return;
843 default:;
844 }
845
846 const bool b_neg = eval_get_sign(x) < 0;
847
848 T xx(x);
849 if (b_neg)
850 xx.negate();
851
852 if (xx.compare(fp_type(0.1)) < 0)
853 {
854 T t1, t2, t3;
855 t1 = ui_type(1);
856 t2 = fp_type(0.5f);
857 t3 = fp_type(1.5f);
858 eval_multiply(xx, xx);
859 xx.negate();
860 hyp2F1(result, t1, t2, t3, xx);
861 eval_multiply(result, x);
862 return;
863 }
864
865 if (xx.compare(fp_type(10)) > 0)
866 {
867 T t1, t2, t3;
868 t1 = fp_type(0.5f);
869 t2 = ui_type(1u);
870 t3 = fp_type(1.5f);
871 eval_multiply(xx, xx);
872 eval_divide(xx, si_type(-1), xx);
873 hyp2F1(result, t1, t2, t3, xx);
874 eval_divide(result, x);
875 if (!b_neg)
876 result.negate();
877 eval_ldexp(t1, get_constant_pi<T>(), -1);
878 eval_add(result, t1);
879 if (b_neg)
880 result.negate();
881 return;
882 }
883
884 // Get initial estimate using standard math function atan.
885 fp_type d;
886 eval_convert_to(&d, xx);
887 result = fp_type(std::atan(d));
888
889 // Newton-Raphson iteration, we should double our precision with each iteration,
890 // in practice this seems to not quite work in all cases... so terminate when we
891 // have at least 2/3 of the digits correct on the assumption that the correction
892 // we've just added will finish the job...
893
894 boost::intmax_t current_precision = eval_ilogb(result);
895 boost::intmax_t target_precision = std::numeric_limits<number<T> >::is_specialized ?
896 current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3
897 : current_precision - 1 - (boost::multiprecision::detail::digits2<number<T> >::value() * 2) / 3;
898
899 T s, c, t;
900 while (current_precision > target_precision)
901 {
902 eval_sin(s, result);
903 eval_cos(c, result);
904 eval_multiply(t, xx, c);
905 eval_subtract(t, s);
906 eval_multiply(s, t, c);
907 eval_add(result, s);
908 current_precision = eval_ilogb(s);
909 if (current_precision <= (std::numeric_limits<typename T::exponent_type>::min)() + 1)
910 break;
911 }
912 if (b_neg)
913 result.negate();
914 }
915
916 template <class T>
eval_atan2(T & result,const T & y,const T & x)917 void eval_atan2(T& result, const T& y, const T& x)
918 {
919 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan2 function is only valid for floating point types.");
920 if (&result == &y)
921 {
922 T temp(y);
923 eval_atan2(result, temp, x);
924 return;
925 }
926 else if (&result == &x)
927 {
928 T temp(x);
929 eval_atan2(result, y, temp);
930 return;
931 }
932
933 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
934
935 switch (eval_fpclassify(y))
936 {
937 case FP_NAN:
938 result = y;
939 errno = EDOM;
940 return;
941 case FP_ZERO:
942 {
943 if (eval_signbit(x))
944 {
945 result = get_constant_pi<T>();
946 if (eval_signbit(y))
947 result.negate();
948 }
949 else
950 {
951 result = y; // Note we allow atan2(0,0) to be +-zero, even though it's mathematically undefined
952 }
953 return;
954 }
955 case FP_INFINITE:
956 {
957 if (eval_fpclassify(x) == FP_INFINITE)
958 {
959 if (eval_signbit(x))
960 {
961 // 3Pi/4
962 eval_ldexp(result, get_constant_pi<T>(), -2);
963 eval_subtract(result, get_constant_pi<T>());
964 if (eval_get_sign(y) >= 0)
965 result.negate();
966 }
967 else
968 {
969 // Pi/4
970 eval_ldexp(result, get_constant_pi<T>(), -2);
971 if (eval_get_sign(y) < 0)
972 result.negate();
973 }
974 }
975 else
976 {
977 eval_ldexp(result, get_constant_pi<T>(), -1);
978 if (eval_get_sign(y) < 0)
979 result.negate();
980 }
981 return;
982 }
983 }
984
985 switch (eval_fpclassify(x))
986 {
987 case FP_NAN:
988 result = x;
989 errno = EDOM;
990 return;
991 case FP_ZERO:
992 {
993 eval_ldexp(result, get_constant_pi<T>(), -1);
994 if (eval_get_sign(y) < 0)
995 result.negate();
996 return;
997 }
998 case FP_INFINITE:
999 if (eval_get_sign(x) > 0)
1000 result = ui_type(0);
1001 else
1002 result = get_constant_pi<T>();
1003 if (eval_get_sign(y) < 0)
1004 result.negate();
1005 return;
1006 }
1007
1008 T xx;
1009 eval_divide(xx, y, x);
1010 if (eval_get_sign(xx) < 0)
1011 xx.negate();
1012
1013 eval_atan(result, xx);
1014
1015 // Determine quadrant (sign) based on signs of x, y
1016 const bool y_neg = eval_get_sign(y) < 0;
1017 const bool x_neg = eval_get_sign(x) < 0;
1018
1019 if (y_neg != x_neg)
1020 result.negate();
1021
1022 if (x_neg)
1023 {
1024 if (y_neg)
1025 eval_subtract(result, get_constant_pi<T>());
1026 else
1027 eval_add(result, get_constant_pi<T>());
1028 }
1029 }
1030 template <class T, class A>
eval_atan2(T & result,const T & x,const A & a)1031 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const T& x, const A& a)
1032 {
1033 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1034 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1035 cast_type c;
1036 c = a;
1037 eval_atan2(result, x, c);
1038 }
1039
1040 template <class T, class A>
eval_atan2(T & result,const A & x,const T & a)1041 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const A& x, const T& a)
1042 {
1043 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1044 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1045 cast_type c;
1046 c = x;
1047 eval_atan2(result, c, a);
1048 }
1049
1050 #ifdef BOOST_MSVC
1051 #pragma warning(pop)
1052 #endif
1053