Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | Download | RSS feed

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