e_acosh.c 1.7 KB

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