annotate DEPENDENCIES/generic/include/boost/multiprecision/detail/functions/trig.hpp @ 133:4acb5d8d80b6 tip

Don't fail environmental check if README.md exists (but .txt and no-suffix don't)
author Chris Cannam
date Tue, 30 Jul 2019 12:25:44 +0100
parents c530137014c0
children
rev   line source
Chris@16 1
Chris@16 2 // Copyright Christopher Kormanyos 2002 - 2011.
Chris@16 3 // Copyright 2011 John Maddock. Distributed under the Boost
Chris@16 4 // Distributed under the Boost Software License, Version 1.0.
Chris@16 5 // (See accompanying file LICENSE_1_0.txt or copy at
Chris@16 6 // http://www.boost.org/LICENSE_1_0.txt)
Chris@16 7
Chris@16 8 // This work is based on an earlier work:
Chris@16 9 // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
Chris@16 10 // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
Chris@16 11 //
Chris@16 12 // This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
Chris@16 13 //
Chris@16 14
Chris@16 15 template <class T>
Chris@16 16 void hyp0F1(T& result, const T& b, const T& x)
Chris@16 17 {
Chris@16 18 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
Chris@16 19 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
Chris@16 20
Chris@16 21 // Compute the series representation of Hypergeometric0F1 taken from
Chris@16 22 // http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F1/06/01/01/
Chris@16 23 // There are no checks on input range or parameter boundaries.
Chris@16 24
Chris@16 25 T x_pow_n_div_n_fact(x);
Chris@16 26 T pochham_b (b);
Chris@16 27 T bp (b);
Chris@16 28
Chris@16 29 eval_divide(result, x_pow_n_div_n_fact, pochham_b);
Chris@16 30 eval_add(result, ui_type(1));
Chris@16 31
Chris@16 32 si_type n;
Chris@16 33
Chris@16 34 T tol;
Chris@16 35 tol = ui_type(1);
Chris@16 36 eval_ldexp(tol, tol, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value);
Chris@16 37 eval_multiply(tol, result);
Chris@16 38 if(eval_get_sign(tol) < 0)
Chris@16 39 tol.negate();
Chris@16 40 T term;
Chris@16 41
Chris@101 42 static const int series_limit =
Chris@16 43 boost::multiprecision::detail::digits2<number<T, et_on> >::value < 100
Chris@16 44 ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value;
Chris@16 45 // Series expansion of hyperg_0f1(; b; x).
Chris@16 46 for(n = 2; n < series_limit; ++n)
Chris@16 47 {
Chris@16 48 eval_multiply(x_pow_n_div_n_fact, x);
Chris@16 49 eval_divide(x_pow_n_div_n_fact, n);
Chris@16 50 eval_increment(bp);
Chris@16 51 eval_multiply(pochham_b, bp);
Chris@16 52
Chris@16 53 eval_divide(term, x_pow_n_div_n_fact, pochham_b);
Chris@16 54 eval_add(result, term);
Chris@16 55
Chris@16 56 bool neg_term = eval_get_sign(term) < 0;
Chris@16 57 if(neg_term)
Chris@16 58 term.negate();
Chris@16 59 if(term.compare(tol) <= 0)
Chris@16 60 break;
Chris@16 61 }
Chris@16 62
Chris@16 63 if(n >= series_limit)
Chris@16 64 BOOST_THROW_EXCEPTION(std::runtime_error("H0F1 Failed to Converge"));
Chris@16 65 }
Chris@16 66
Chris@16 67
Chris@16 68 template <class T>
Chris@16 69 void eval_sin(T& result, const T& x)
Chris@16 70 {
Chris@16 71 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sin function is only valid for floating point types.");
Chris@16 72 if(&result == &x)
Chris@16 73 {
Chris@16 74 T temp;
Chris@16 75 eval_sin(temp, x);
Chris@16 76 result = temp;
Chris@16 77 return;
Chris@16 78 }
Chris@16 79
Chris@16 80 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
Chris@16 81 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
Chris@16 82 typedef typename mpl::front<typename T::float_types>::type fp_type;
Chris@16 83
Chris@16 84 switch(eval_fpclassify(x))
Chris@16 85 {
Chris@16 86 case FP_INFINITE:
Chris@16 87 case FP_NAN:
Chris@16 88 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
Chris@16 89 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
Chris@16 90 else
Chris@16 91 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
Chris@16 92 return;
Chris@16 93 case FP_ZERO:
Chris@16 94 result = ui_type(0);
Chris@16 95 return;
Chris@16 96 default: ;
Chris@16 97 }
Chris@16 98
Chris@16 99 // Local copy of the argument
Chris@16 100 T xx = x;
Chris@16 101
Chris@16 102 // Analyze and prepare the phase of the argument.
Chris@16 103 // Make a local, positive copy of the argument, xx.
Chris@16 104 // The argument xx will be reduced to 0 <= xx <= pi/2.
Chris@16 105 bool b_negate_sin = false;
Chris@16 106
Chris@16 107 if(eval_get_sign(x) < 0)
Chris@16 108 {
Chris@16 109 xx.negate();
Chris@16 110 b_negate_sin = !b_negate_sin;
Chris@16 111 }
Chris@16 112
Chris@16 113 T n_pi, t;
Chris@16 114 // Remove even multiples of pi.
Chris@16 115 if(xx.compare(get_constant_pi<T>()) > 0)
Chris@16 116 {
Chris@16 117 eval_divide(n_pi, xx, get_constant_pi<T>());
Chris@16 118 eval_trunc(n_pi, n_pi);
Chris@16 119 t = ui_type(2);
Chris@16 120 eval_fmod(t, n_pi, t);
Chris@16 121 const bool b_n_pi_is_even = eval_get_sign(t) == 0;
Chris@16 122 eval_multiply(n_pi, get_constant_pi<T>());
Chris@16 123 eval_subtract(xx, n_pi);
Chris@16 124
Chris@16 125 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
Chris@16 126 BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
Chris@16 127
Chris@16 128 // Adjust signs if the multiple of pi is not even.
Chris@16 129 if(!b_n_pi_is_even)
Chris@16 130 {
Chris@16 131 b_negate_sin = !b_negate_sin;
Chris@16 132 }
Chris@16 133 }
Chris@16 134
Chris@16 135 // Reduce the argument to 0 <= xx <= pi/2.
Chris@16 136 eval_ldexp(t, get_constant_pi<T>(), -1);
Chris@16 137 if(xx.compare(t) > 0)
Chris@16 138 {
Chris@16 139 eval_subtract(xx, get_constant_pi<T>(), xx);
Chris@16 140 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
Chris@16 141 }
Chris@16 142
Chris@16 143 eval_subtract(t, xx);
Chris@16 144 const bool b_zero = eval_get_sign(xx) == 0;
Chris@16 145 const bool b_pi_half = eval_get_sign(t) == 0;
Chris@16 146
Chris@16 147 // Check if the reduced argument is very close to 0 or pi/2.
Chris@16 148 const bool b_near_zero = xx.compare(fp_type(1e-1)) < 0;
Chris@16 149 const bool b_near_pi_half = t.compare(fp_type(1e-1)) < 0;;
Chris@16 150
Chris@16 151 if(b_zero)
Chris@16 152 {
Chris@16 153 result = ui_type(0);
Chris@16 154 }
Chris@16 155 else if(b_pi_half)
Chris@16 156 {
Chris@16 157 result = ui_type(1);
Chris@16 158 }
Chris@16 159 else if(b_near_zero)
Chris@16 160 {
Chris@16 161 eval_multiply(t, xx, xx);
Chris@16 162 eval_divide(t, si_type(-4));
Chris@16 163 T t2;
Chris@16 164 t2 = fp_type(1.5);
Chris@16 165 hyp0F1(result, t2, t);
Chris@16 166 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
Chris@16 167 eval_multiply(result, xx);
Chris@16 168 }
Chris@16 169 else if(b_near_pi_half)
Chris@16 170 {
Chris@16 171 eval_multiply(t, t);
Chris@16 172 eval_divide(t, si_type(-4));
Chris@16 173 T t2;
Chris@16 174 t2 = fp_type(0.5);
Chris@16 175 hyp0F1(result, t2, t);
Chris@16 176 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
Chris@16 177 }
Chris@16 178 else
Chris@16 179 {
Chris@16 180 // Scale to a small argument for an efficient Taylor series,
Chris@16 181 // implemented as a hypergeometric function. Use a standard
Chris@16 182 // divide by three identity a certain number of times.
Chris@16 183 // Here we use division by 3^9 --> (19683 = 3^9).
Chris@16 184
Chris@16 185 static const si_type n_scale = 9;
Chris@16 186 static const si_type n_three_pow_scale = static_cast<si_type>(19683L);
Chris@16 187
Chris@16 188 eval_divide(xx, n_three_pow_scale);
Chris@16 189
Chris@16 190 // Now with small arguments, we are ready for a series expansion.
Chris@16 191 eval_multiply(t, xx, xx);
Chris@16 192 eval_divide(t, si_type(-4));
Chris@16 193 T t2;
Chris@16 194 t2 = fp_type(1.5);
Chris@16 195 hyp0F1(result, t2, t);
Chris@16 196 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
Chris@16 197 eval_multiply(result, xx);
Chris@16 198
Chris@16 199 // Convert back using multiple angle identity.
Chris@16 200 for(boost::int32_t k = static_cast<boost::int32_t>(0); k < n_scale; k++)
Chris@16 201 {
Chris@16 202 // Rescale the cosine value using the multiple angle identity.
Chris@16 203 eval_multiply(t2, result, ui_type(3));
Chris@16 204 eval_multiply(t, result, result);
Chris@16 205 eval_multiply(t, result);
Chris@16 206 eval_multiply(t, ui_type(4));
Chris@16 207 eval_subtract(result, t2, t);
Chris@16 208 }
Chris@16 209 }
Chris@16 210
Chris@16 211 if(b_negate_sin)
Chris@16 212 result.negate();
Chris@16 213 }
Chris@16 214
Chris@16 215 template <class T>
Chris@16 216 void eval_cos(T& result, const T& x)
Chris@16 217 {
Chris@16 218 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cos function is only valid for floating point types.");
Chris@16 219 if(&result == &x)
Chris@16 220 {
Chris@16 221 T temp;
Chris@16 222 eval_cos(temp, x);
Chris@16 223 result = temp;
Chris@16 224 return;
Chris@16 225 }
Chris@16 226
Chris@16 227 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
Chris@16 228 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
Chris@16 229 typedef typename mpl::front<typename T::float_types>::type fp_type;
Chris@16 230
Chris@16 231 switch(eval_fpclassify(x))
Chris@16 232 {
Chris@16 233 case FP_INFINITE:
Chris@16 234 case FP_NAN:
Chris@16 235 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
Chris@16 236 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
Chris@16 237 else
Chris@16 238 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
Chris@16 239 return;
Chris@16 240 case FP_ZERO:
Chris@16 241 result = ui_type(1);
Chris@16 242 return;
Chris@16 243 default: ;
Chris@16 244 }
Chris@16 245
Chris@16 246 // Local copy of the argument
Chris@16 247 T xx = x;
Chris@16 248
Chris@16 249 // Analyze and prepare the phase of the argument.
Chris@16 250 // Make a local, positive copy of the argument, xx.
Chris@16 251 // The argument xx will be reduced to 0 <= xx <= pi/2.
Chris@16 252 bool b_negate_cos = false;
Chris@16 253
Chris@16 254 if(eval_get_sign(x) < 0)
Chris@16 255 {
Chris@16 256 xx.negate();
Chris@16 257 }
Chris@16 258
Chris@16 259 T n_pi, t;
Chris@16 260 // Remove even multiples of pi.
Chris@16 261 if(xx.compare(get_constant_pi<T>()) > 0)
Chris@16 262 {
Chris@16 263 eval_divide(t, xx, get_constant_pi<T>());
Chris@16 264 eval_trunc(n_pi, t);
Chris@16 265 BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
Chris@16 266 eval_multiply(t, n_pi, get_constant_pi<T>());
Chris@16 267 BOOST_MATH_INSTRUMENT_CODE(t.str(0, std::ios_base::scientific));
Chris@16 268 eval_subtract(xx, t);
Chris@16 269 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
Chris@16 270
Chris@16 271 // Adjust signs if the multiple of pi is not even.
Chris@16 272 t = ui_type(2);
Chris@16 273 eval_fmod(t, n_pi, t);
Chris@16 274 const bool b_n_pi_is_even = eval_get_sign(t) == 0;
Chris@16 275
Chris@16 276 if(!b_n_pi_is_even)
Chris@16 277 {
Chris@16 278 b_negate_cos = !b_negate_cos;
Chris@16 279 }
Chris@16 280 }
Chris@16 281
Chris@16 282 // Reduce the argument to 0 <= xx <= pi/2.
Chris@16 283 eval_ldexp(t, get_constant_pi<T>(), -1);
Chris@16 284 int com = xx.compare(t);
Chris@16 285 if(com > 0)
Chris@16 286 {
Chris@16 287 eval_subtract(xx, get_constant_pi<T>(), xx);
Chris@16 288 b_negate_cos = !b_negate_cos;
Chris@16 289 BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
Chris@16 290 }
Chris@16 291
Chris@16 292 const bool b_zero = eval_get_sign(xx) == 0;
Chris@16 293 const bool b_pi_half = com == 0;
Chris@16 294
Chris@16 295 // Check if the reduced argument is very close to 0.
Chris@16 296 const bool b_near_zero = xx.compare(fp_type(1e-1)) < 0;
Chris@16 297
Chris@16 298 if(b_zero)
Chris@16 299 {
Chris@16 300 result = si_type(1);
Chris@16 301 }
Chris@16 302 else if(b_pi_half)
Chris@16 303 {
Chris@16 304 result = si_type(0);
Chris@16 305 }
Chris@16 306 else if(b_near_zero)
Chris@16 307 {
Chris@16 308 eval_multiply(t, xx, xx);
Chris@16 309 eval_divide(t, si_type(-4));
Chris@16 310 n_pi = fp_type(0.5f);
Chris@16 311 hyp0F1(result, n_pi, t);
Chris@16 312 BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
Chris@16 313 }
Chris@16 314 else
Chris@16 315 {
Chris@16 316 eval_subtract(t, xx);
Chris@16 317 eval_sin(result, t);
Chris@16 318 }
Chris@16 319 if(b_negate_cos)
Chris@16 320 result.negate();
Chris@16 321 }
Chris@16 322
Chris@16 323 template <class T>
Chris@16 324 void eval_tan(T& result, const T& x)
Chris@16 325 {
Chris@16 326 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tan function is only valid for floating point types.");
Chris@16 327 if(&result == &x)
Chris@16 328 {
Chris@16 329 T temp;
Chris@16 330 eval_tan(temp, x);
Chris@16 331 result = temp;
Chris@16 332 return;
Chris@16 333 }
Chris@16 334 T t;
Chris@16 335 eval_sin(result, x);
Chris@16 336 eval_cos(t, x);
Chris@16 337 eval_divide(result, t);
Chris@16 338 }
Chris@16 339
Chris@16 340 template <class T>
Chris@16 341 void hyp2F1(T& result, const T& a, const T& b, const T& c, const T& x)
Chris@16 342 {
Chris@16 343 // Compute the series representation of hyperg_2f1 taken from
Chris@16 344 // Abramowitz and Stegun 15.1.1.
Chris@16 345 // There are no checks on input range or parameter boundaries.
Chris@16 346
Chris@16 347 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
Chris@16 348
Chris@16 349 T x_pow_n_div_n_fact(x);
Chris@16 350 T pochham_a (a);
Chris@16 351 T pochham_b (b);
Chris@16 352 T pochham_c (c);
Chris@16 353 T ap (a);
Chris@16 354 T bp (b);
Chris@16 355 T cp (c);
Chris@16 356
Chris@16 357 eval_multiply(result, pochham_a, pochham_b);
Chris@16 358 eval_divide(result, pochham_c);
Chris@16 359 eval_multiply(result, x_pow_n_div_n_fact);
Chris@16 360 eval_add(result, ui_type(1));
Chris@16 361
Chris@16 362 T lim;
Chris@16 363 eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value);
Chris@16 364
Chris@16 365 if(eval_get_sign(lim) < 0)
Chris@16 366 lim.negate();
Chris@16 367
Chris@16 368 ui_type n;
Chris@16 369 T term;
Chris@16 370
Chris@16 371 static const unsigned series_limit =
Chris@16 372 boost::multiprecision::detail::digits2<number<T, et_on> >::value < 100
Chris@16 373 ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value;
Chris@16 374 // Series expansion of hyperg_2f1(a, b; c; x).
Chris@16 375 for(n = 2; n < series_limit; ++n)
Chris@16 376 {
Chris@16 377 eval_multiply(x_pow_n_div_n_fact, x);
Chris@16 378 eval_divide(x_pow_n_div_n_fact, n);
Chris@16 379
Chris@16 380 eval_increment(ap);
Chris@16 381 eval_multiply(pochham_a, ap);
Chris@16 382 eval_increment(bp);
Chris@16 383 eval_multiply(pochham_b, bp);
Chris@16 384 eval_increment(cp);
Chris@16 385 eval_multiply(pochham_c, cp);
Chris@16 386
Chris@16 387 eval_multiply(term, pochham_a, pochham_b);
Chris@16 388 eval_divide(term, pochham_c);
Chris@16 389 eval_multiply(term, x_pow_n_div_n_fact);
Chris@16 390 eval_add(result, term);
Chris@16 391
Chris@16 392 if(eval_get_sign(term) < 0)
Chris@16 393 term.negate();
Chris@16 394 if(lim.compare(term) >= 0)
Chris@16 395 break;
Chris@16 396 }
Chris@16 397 if(n > series_limit)
Chris@16 398 BOOST_THROW_EXCEPTION(std::runtime_error("H2F1 failed to converge."));
Chris@16 399 }
Chris@16 400
Chris@16 401 template <class T>
Chris@16 402 void eval_asin(T& result, const T& x)
Chris@16 403 {
Chris@16 404 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The asin function is only valid for floating point types.");
Chris@16 405 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
Chris@16 406 typedef typename mpl::front<typename T::float_types>::type fp_type;
Chris@16 407
Chris@16 408 if(&result == &x)
Chris@16 409 {
Chris@16 410 T t(x);
Chris@16 411 eval_asin(result, t);
Chris@16 412 return;
Chris@16 413 }
Chris@16 414
Chris@16 415 switch(eval_fpclassify(x))
Chris@16 416 {
Chris@16 417 case FP_NAN:
Chris@16 418 case FP_INFINITE:
Chris@16 419 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
Chris@16 420 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
Chris@16 421 else
Chris@16 422 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
Chris@16 423 return;
Chris@16 424 case FP_ZERO:
Chris@16 425 result = ui_type(0);
Chris@16 426 return;
Chris@16 427 default: ;
Chris@16 428 }
Chris@16 429
Chris@16 430 const bool b_neg = eval_get_sign(x) < 0;
Chris@16 431
Chris@16 432 T xx(x);
Chris@16 433 if(b_neg)
Chris@16 434 xx.negate();
Chris@16 435
Chris@16 436 int c = xx.compare(ui_type(1));
Chris@16 437 if(c > 0)
Chris@16 438 {
Chris@16 439 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
Chris@16 440 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
Chris@16 441 else
Chris@16 442 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
Chris@16 443 return;
Chris@16 444 }
Chris@16 445 else if(c == 0)
Chris@16 446 {
Chris@16 447 result = get_constant_pi<T>();
Chris@16 448 eval_ldexp(result, result, -1);
Chris@16 449 if(b_neg)
Chris@16 450 result.negate();
Chris@16 451 return;
Chris@16 452 }
Chris@16 453
Chris@16 454 if(xx.compare(fp_type(1e-4)) < 0)
Chris@16 455 {
Chris@16 456 // http://functions.wolfram.com/ElementaryFunctions/ArcSin/26/01/01/
Chris@16 457 eval_multiply(xx, xx);
Chris@16 458 T t1, t2;
Chris@16 459 t1 = fp_type(0.5f);
Chris@16 460 t2 = fp_type(1.5f);
Chris@16 461 hyp2F1(result, t1, t1, t2, xx);
Chris@16 462 eval_multiply(result, x);
Chris@16 463 return;
Chris@16 464 }
Chris@16 465 else if(xx.compare(fp_type(1 - 1e-4f)) > 0)
Chris@16 466 {
Chris@16 467 T dx1;
Chris@16 468 T t1, t2;
Chris@16 469 eval_subtract(dx1, ui_type(1), xx);
Chris@16 470 t1 = fp_type(0.5f);
Chris@16 471 t2 = fp_type(1.5f);
Chris@16 472 eval_ldexp(dx1, dx1, -1);
Chris@16 473 hyp2F1(result, t1, t1, t2, dx1);
Chris@16 474 eval_ldexp(dx1, dx1, 2);
Chris@16 475 eval_sqrt(t1, dx1);
Chris@16 476 eval_multiply(result, t1);
Chris@16 477 eval_ldexp(t1, get_constant_pi<T>(), -1);
Chris@16 478 result.negate();
Chris@16 479 eval_add(result, t1);
Chris@16 480 if(b_neg)
Chris@16 481 result.negate();
Chris@16 482 return;
Chris@16 483 }
Chris@16 484
Chris@16 485 // Get initial estimate using standard math function asin.
Chris@16 486 double dd;
Chris@16 487 eval_convert_to(&dd, xx);
Chris@16 488
Chris@16 489 result = fp_type(std::asin(dd));
Chris@16 490
Chris@101 491 unsigned current_digits = std::numeric_limits<double>::digits - 5;
Chris@101 492 unsigned target_precision = boost::multiprecision::detail::digits2<number<T, et_on> >::value;
Chris@101 493
Chris@16 494 // Newton-Raphson iteration
Chris@101 495 while(current_digits < target_precision)
Chris@16 496 {
Chris@16 497 T s, c;
Chris@16 498 eval_sin(s, result);
Chris@16 499 eval_cos(c, result);
Chris@16 500 eval_subtract(s, xx);
Chris@16 501 eval_divide(s, c);
Chris@16 502 eval_subtract(result, s);
Chris@16 503
Chris@101 504 current_digits *= 2;
Chris@101 505 /*
Chris@16 506 T lim;
Chris@16 507 eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value);
Chris@16 508 if(eval_get_sign(s) < 0)
Chris@16 509 s.negate();
Chris@16 510 if(eval_get_sign(lim) < 0)
Chris@16 511 lim.negate();
Chris@16 512 if(lim.compare(s) >= 0)
Chris@16 513 break;
Chris@101 514 */
Chris@16 515 }
Chris@16 516 if(b_neg)
Chris@16 517 result.negate();
Chris@16 518 }
Chris@16 519
Chris@16 520 template <class T>
Chris@16 521 inline void eval_acos(T& result, const T& x)
Chris@16 522 {
Chris@16 523 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The acos function is only valid for floating point types.");
Chris@16 524 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
Chris@16 525
Chris@16 526 switch(eval_fpclassify(x))
Chris@16 527 {
Chris@16 528 case FP_NAN:
Chris@16 529 case FP_INFINITE:
Chris@16 530 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
Chris@16 531 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
Chris@16 532 else
Chris@16 533 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
Chris@16 534 return;
Chris@16 535 case FP_ZERO:
Chris@16 536 result = get_constant_pi<T>();
Chris@16 537 eval_ldexp(result, result, -1); // divide by two.
Chris@16 538 return;
Chris@16 539 }
Chris@16 540
Chris@16 541 eval_abs(result, x);
Chris@16 542 int c = result.compare(ui_type(1));
Chris@16 543
Chris@16 544 if(c > 0)
Chris@16 545 {
Chris@16 546 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
Chris@16 547 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
Chris@16 548 else
Chris@16 549 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
Chris@16 550 return;
Chris@16 551 }
Chris@16 552 else if(c == 0)
Chris@16 553 {
Chris@16 554 if(eval_get_sign(x) < 0)
Chris@16 555 result = get_constant_pi<T>();
Chris@16 556 else
Chris@16 557 result = ui_type(0);
Chris@16 558 return;
Chris@16 559 }
Chris@16 560
Chris@16 561 eval_asin(result, x);
Chris@16 562 T t;
Chris@16 563 eval_ldexp(t, get_constant_pi<T>(), -1);
Chris@16 564 eval_subtract(result, t);
Chris@16 565 result.negate();
Chris@16 566 }
Chris@16 567
Chris@16 568 template <class T>
Chris@16 569 void eval_atan(T& result, const T& x)
Chris@16 570 {
Chris@16 571 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan function is only valid for floating point types.");
Chris@16 572 typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
Chris@16 573 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
Chris@16 574 typedef typename mpl::front<typename T::float_types>::type fp_type;
Chris@16 575
Chris@16 576 switch(eval_fpclassify(x))
Chris@16 577 {
Chris@16 578 case FP_NAN:
Chris@16 579 result = x;
Chris@16 580 return;
Chris@16 581 case FP_ZERO:
Chris@16 582 result = ui_type(0);
Chris@16 583 return;
Chris@16 584 case FP_INFINITE:
Chris@16 585 if(eval_get_sign(x) < 0)
Chris@16 586 {
Chris@16 587 eval_ldexp(result, get_constant_pi<T>(), -1);
Chris@16 588 result.negate();
Chris@16 589 }
Chris@16 590 else
Chris@16 591 eval_ldexp(result, get_constant_pi<T>(), -1);
Chris@16 592 return;
Chris@16 593 default: ;
Chris@16 594 }
Chris@16 595
Chris@16 596 const bool b_neg = eval_get_sign(x) < 0;
Chris@16 597
Chris@16 598 T xx(x);
Chris@16 599 if(b_neg)
Chris@16 600 xx.negate();
Chris@16 601
Chris@16 602 if(xx.compare(fp_type(0.1)) < 0)
Chris@16 603 {
Chris@16 604 T t1, t2, t3;
Chris@16 605 t1 = ui_type(1);
Chris@16 606 t2 = fp_type(0.5f);
Chris@16 607 t3 = fp_type(1.5f);
Chris@16 608 eval_multiply(xx, xx);
Chris@16 609 xx.negate();
Chris@16 610 hyp2F1(result, t1, t2, t3, xx);
Chris@16 611 eval_multiply(result, x);
Chris@16 612 return;
Chris@16 613 }
Chris@16 614
Chris@16 615 if(xx.compare(fp_type(10)) > 0)
Chris@16 616 {
Chris@16 617 T t1, t2, t3;
Chris@16 618 t1 = fp_type(0.5f);
Chris@16 619 t2 = ui_type(1u);
Chris@16 620 t3 = fp_type(1.5f);
Chris@16 621 eval_multiply(xx, xx);
Chris@16 622 eval_divide(xx, si_type(-1), xx);
Chris@16 623 hyp2F1(result, t1, t2, t3, xx);
Chris@16 624 eval_divide(result, x);
Chris@16 625 if(!b_neg)
Chris@16 626 result.negate();
Chris@16 627 eval_ldexp(t1, get_constant_pi<T>(), -1);
Chris@16 628 eval_add(result, t1);
Chris@16 629 if(b_neg)
Chris@16 630 result.negate();
Chris@16 631 return;
Chris@16 632 }
Chris@16 633
Chris@16 634
Chris@16 635 // Get initial estimate using standard math function atan.
Chris@16 636 fp_type d;
Chris@16 637 eval_convert_to(&d, xx);
Chris@16 638 result = fp_type(std::atan(d));
Chris@16 639
Chris@16 640 // Newton-Raphson iteration
Chris@16 641 static const boost::int32_t double_digits10_minus_a_few = std::numeric_limits<double>::digits10 - 3;
Chris@16 642
Chris@16 643 T s, c, t;
Chris@16 644 for(boost::int32_t digits = double_digits10_minus_a_few; digits <= std::numeric_limits<number<T, et_on> >::digits10; digits *= 2)
Chris@16 645 {
Chris@16 646 eval_sin(s, result);
Chris@16 647 eval_cos(c, result);
Chris@16 648 eval_multiply(t, xx, c);
Chris@16 649 eval_subtract(t, s);
Chris@16 650 eval_multiply(s, t, c);
Chris@16 651 eval_add(result, s);
Chris@16 652 }
Chris@16 653 if(b_neg)
Chris@16 654 result.negate();
Chris@16 655 }
Chris@16 656
Chris@16 657 template <class T>
Chris@16 658 void eval_atan2(T& result, const T& y, const T& x)
Chris@16 659 {
Chris@16 660 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan2 function is only valid for floating point types.");
Chris@16 661 if(&result == &y)
Chris@16 662 {
Chris@16 663 T temp(y);
Chris@16 664 eval_atan2(result, temp, x);
Chris@16 665 return;
Chris@16 666 }
Chris@16 667 else if(&result == &x)
Chris@16 668 {
Chris@16 669 T temp(x);
Chris@16 670 eval_atan2(result, y, temp);
Chris@16 671 return;
Chris@16 672 }
Chris@16 673
Chris@16 674 typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
Chris@16 675
Chris@16 676 switch(eval_fpclassify(y))
Chris@16 677 {
Chris@16 678 case FP_NAN:
Chris@16 679 result = y;
Chris@16 680 return;
Chris@16 681 case FP_ZERO:
Chris@16 682 {
Chris@16 683 int c = eval_get_sign(x);
Chris@16 684 if(c < 0)
Chris@16 685 result = get_constant_pi<T>();
Chris@16 686 else if(c >= 0)
Chris@16 687 result = ui_type(0); // Note we allow atan2(0,0) to be zero, even though it's mathematically undefined
Chris@16 688 return;
Chris@16 689 }
Chris@16 690 case FP_INFINITE:
Chris@16 691 {
Chris@16 692 if(eval_fpclassify(x) == FP_INFINITE)
Chris@16 693 {
Chris@16 694 if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
Chris@16 695 result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
Chris@16 696 else
Chris@16 697 BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
Chris@16 698 }
Chris@16 699 else
Chris@16 700 {
Chris@16 701 eval_ldexp(result, get_constant_pi<T>(), -1);
Chris@16 702 if(eval_get_sign(y) < 0)
Chris@16 703 result.negate();
Chris@16 704 }
Chris@16 705 return;
Chris@16 706 }
Chris@16 707 }
Chris@16 708
Chris@16 709 switch(eval_fpclassify(x))
Chris@16 710 {
Chris@16 711 case FP_NAN:
Chris@16 712 result = x;
Chris@16 713 return;
Chris@16 714 case FP_ZERO:
Chris@16 715 {
Chris@16 716 eval_ldexp(result, get_constant_pi<T>(), -1);
Chris@16 717 if(eval_get_sign(y) < 0)
Chris@16 718 result.negate();
Chris@16 719 return;
Chris@16 720 }
Chris@16 721 case FP_INFINITE:
Chris@16 722 if(eval_get_sign(x) > 0)
Chris@16 723 result = ui_type(0);
Chris@16 724 else
Chris@16 725 result = get_constant_pi<T>();
Chris@16 726 if(eval_get_sign(y) < 0)
Chris@16 727 result.negate();
Chris@16 728 return;
Chris@16 729 }
Chris@16 730
Chris@16 731 T xx;
Chris@16 732 eval_divide(xx, y, x);
Chris@16 733 if(eval_get_sign(xx) < 0)
Chris@16 734 xx.negate();
Chris@16 735
Chris@16 736 eval_atan(result, xx);
Chris@16 737
Chris@16 738 // Determine quadrant (sign) based on signs of x, y
Chris@16 739 const bool y_neg = eval_get_sign(y) < 0;
Chris@16 740 const bool x_neg = eval_get_sign(x) < 0;
Chris@16 741
Chris@16 742 if(y_neg != x_neg)
Chris@16 743 result.negate();
Chris@16 744
Chris@16 745 if(x_neg)
Chris@16 746 {
Chris@16 747 if(y_neg)
Chris@16 748 eval_subtract(result, get_constant_pi<T>());
Chris@16 749 else
Chris@16 750 eval_add(result, get_constant_pi<T>());
Chris@16 751 }
Chris@16 752 }
Chris@16 753 template<class T, class A>
Chris@16 754 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const T& x, const A& a)
Chris@16 755 {
Chris@16 756 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
Chris@16 757 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
Chris@16 758 cast_type c;
Chris@16 759 c = a;
Chris@16 760 eval_atan2(result, x, c);
Chris@16 761 }
Chris@16 762
Chris@16 763 template<class T, class A>
Chris@16 764 inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const A& x, const T& a)
Chris@16 765 {
Chris@16 766 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
Chris@16 767 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
Chris@16 768 cast_type c;
Chris@16 769 c = x;
Chris@16 770 eval_atan2(result, c, a);
Chris@16 771 }
Chris@16 772