sysdeps/ieee754/ldbl-96/e_acoshl.c

   1 /* e_acoshl.c -- long double version of e_acosh.c.
   2  * Conversion to long double by Ulrich Drepper,
   3  * Cygnus Support, drepper@cygnus.com.
   4  */
   5
   6 /*
   7  * ====================================================
   8  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   9  *
  10  * Developed at SunPro, a Sun Microsystems, Inc. business.
  11  * Permission to use, copy, modify, and distribute this
  12  * software is freely granted, provided that this notice
  13  * is preserved.
  14  * ====================================================
  15  */
  16
  17 /* __ieee754_acoshl(x)
  18  * Method :
  19  *      Based on
  20  *              acoshl(x) = logl [ x + sqrtl(x*x-1) ]
  21  *      we have
  22  *              acoshl(x) := logl(x)+ln2,       if x is large; else
  23  *              acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
  24  *              acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
  25  *
  26  * Special cases:
  27  *      acoshl(x) is NaN with signal if x<1.
  28  *      acoshl(NaN) is NaN without signal.
  29  */
  30
  31 #include <math.h>
  32 #include <math_private.h>
  33
  34 static const long double
  35 one     = 1.0,
  36 ln2     = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
  37
  38 long double
  39 __ieee754_acoshl(long double x)
  40 {
  41         long double t;
  42         u_int32_t se,i0,i1;
  43         GET_LDOUBLE_WORDS(se,i0,i1,x);
  44         if(se<0x3fff || se & 0x8000) {  /* x < 1 */
  45             return (x-x)/(x-x);
  46         } else if(se >=0x401d) {        /* x > 2**30 */
  47             if(se >=0x7fff) {           /* x is inf of NaN */
  48                 return x+x;
  49             } else
  50                 return __ieee754_logl(x)+ln2;   /* acoshl(huge)=logl(2x) */
  51         } else if(((se-0x3fff)|(i0^0x80000000)|i1)==0) {
  52             return 0.0;                 /* acosh(1) = 0 */
  53         } else if (se > 0x4000) {       /* 2**28 > x > 2 */
  54             t=x*x;
  55             return __ieee754_logl(2.0*x-one/(x+__ieee754_sqrtl(t-one)));
  56         } else {                        /* 1<x<2 */
  57             t = x-one;
  58             return __log1pl(t+__ieee754_sqrtl(2.0*t+t*t));
  59         }
  60 }
  61 strong_alias (__ieee754_acoshl, __acoshl_finite)