e_acosh.c 1.8 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. /*
  2. * ====================================================
  3. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  4. *
  5. * Developed at SunPro, a Sun Microsystems, Inc. business.
  6. * Permission to use, copy, modify, and distribute this
  7. * software is freely granted, provided that this notice
  8. * is preserved.
  9. * ====================================================
  10. */
  11. /* __ieee754_acosh(x)
  12. * Method :
  13. * Based on
  14. * acosh(x) = log [ x + sqrt(x*x-1) ]
  15. * we have
  16. * acosh(x) := log(x)+ln2, if x is large; else
  17. * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
  18. * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
  19. *
  20. * Special cases:
  21. * acosh(x) is NaN with signal if x<1.
  22. * acosh(NaN) is NaN without signal.
  23. */
  24. #include "math.h"
  25. #include "math_private.h"
  26. static const double
  27. one = 1.0,
  28. ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
  29. double __ieee754_acosh(double x)
  30. {
  31. double t;
  32. int32_t hx;
  33. u_int32_t lx;
  34. EXTRACT_WORDS(hx,lx,x);
  35. if(hx<0x3ff00000) { /* x < 1 */
  36. return (x-x)/(x-x);
  37. } else if(hx >=0x41b00000) { /* x > 2**28 */
  38. if(hx >=0x7ff00000) { /* x is inf of NaN */
  39. return x+x;
  40. } else
  41. return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */
  42. } else if(((hx-0x3ff00000)|lx)==0) {
  43. return 0.0; /* acosh(1) = 0 */
  44. } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
  45. t=x*x;
  46. return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one)));
  47. } else { /* 1<x<2 */
  48. t = x-one;
  49. return log1p(t+sqrt(2.0*t+t*t));
  50. }
  51. }
  52. /*
  53. * wrapper acosh(x)
  54. */
  55. #ifndef _IEEE_LIBM
  56. double acosh(double x)
  57. {
  58. double z = __ieee754_acosh(x);
  59. if (_LIB_VERSION == _IEEE_ || isnan(x))
  60. return z;
  61. if (x < 1.0)
  62. return __kernel_standard(x, x, 29); /* acosh(x<1) */
  63. return z;
  64. }
  65. #else
  66. strong_alias(__ieee754_acosh, acosh)
  67. #endif
  68. libm_hidden_def(acosh)