Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1.  
  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.  
  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. #ifndef _DOUBLE_IS_32BITS
  32.  
  33. #ifdef __STDC__
  34. static const double
  35. #else
  36. static double
  37. #endif
  38. one     = 1.0,
  39. ln2     = 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
  40.  
  41. #ifdef __STDC__
  42.         double __ieee754_acosh(double x)
  43. #else
  44.         double __ieee754_acosh(x)
  45.         double x;
  46. #endif
  47. {      
  48.         double t;
  49.         __int32_t hx;
  50.         __uint32_t lx;
  51.         EXTRACT_WORDS(hx,lx,x);
  52.         if(hx<0x3ff00000) {             /* x < 1 */
  53.             return (x-x)/(x-x);
  54.         } else if(hx >=0x41b00000) {    /* x > 2**28 */
  55.             if(hx >=0x7ff00000) {       /* x is inf of NaN */
  56.                 return x+x;
  57.             } else
  58.                 return __ieee754_log(x)+ln2;    /* acosh(huge)=log(2x) */
  59.         } else if(((hx-0x3ff00000)|lx)==0) {
  60.             return 0.0;                 /* acosh(1) = 0 */
  61.         } else if (hx > 0x40000000) {   /* 2**28 > x > 2 */
  62.             t=x*x;
  63.             return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one)));
  64.         } else {                        /* 1<x<2 */
  65.             t = x-one;
  66.             return log1p(t+__ieee754_sqrt(2.0*t+t*t));
  67.         }
  68. }
  69.  
  70. #endif /* defined(_DOUBLE_IS_32BITS) */
  71.