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

source code of boost/libs/multiprecision/include/boost/multiprecision/detail/functions/trig.hpp