home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / math / e_acosh.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-16  |  1.6 KB  |  66 lines
  1.  
  2. /* @(#)e_acosh.c 1.3 95/01/18 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  6.  *
  7.  * Developed at SunSoft, 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.  
  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.  
  29. #include "fdlibm.h"
  30.  
  31. #ifdef __STDC__
  32. static const double 
  33. #else
  34. static double 
  35. #endif
  36. one    = 1.0,
  37. ln2    = 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
  38.  
  39. #ifdef __STDC__
  40.     double __ieee754_acosh(double x)
  41. #else
  42.     double __ieee754_acosh(x)
  43.     double x;
  44. #endif
  45. {    
  46.     double t;
  47.     int hx;
  48.     hx = __HI(x);
  49.     if(hx<0x3ff00000) {        /* x < 1 */
  50.         return (x-x)/(x-x);
  51.     } else if(hx >=0x41b00000) {    /* x > 2**28 */
  52.         if(hx >=0x7ff00000) {    /* x is inf of NaN */
  53.             return x+x;
  54.         } else 
  55.         return __ieee754_log(x)+ln2;    /* acosh(huge)=log(2x) */
  56.     } else if(((hx-0x3ff00000)|__LO(x))==0) {
  57.         return 0.0;            /* acosh(1) = 0 */
  58.     } else if (hx > 0x40000000) {    /* 2**28 > x > 2 */
  59.         t=x*x;
  60.         return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
  61.     } else {            /* 1<x<2 */
  62.         t = x-one;
  63.         return log1p(t+sqrt(2.0*t+t*t));
  64.     }
  65. }
  66.