sysdeps/libm-ieee754/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 #if defined(LIBM_SCCS) && !defined(lint)
  18 static char rcsid[] = "$NetBSD: $";
  19 #endif
  20
  21 /* __ieee754_acoshl(x)
  22  * Method :
  23  *      Based on
  24  *              acoshl(x) = logl [ x + sqrtl(x*x-1) ]
  25  *      we have
  26  *              acoshl(x) := logl(x)+ln2,       if x is large; else
  27  *              acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
  28  *              acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
  29  *
  30  * Special cases:
  31  *      acoshl(x) is NaN with signal if x<1.
  32  *      acoshl(NaN) is NaN without signal.
  33  */
  34
  35 #include "math.h"
  36 #include "math_private.h"
  37
  38 #ifdef __STDC__
  39 static const long double
  40 #else
  41 static long double
  42 #endif
  43 one     = 1.0,
  44 ln2     = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
  45
  46 #ifdef __STDC__
  47         long double __ieee754_acoshl(long double x)
  48 #else
  49         long double __ieee754_acoshl(x)
  50         long double x;
  51 #endif
  52 {
  53         long double t;
  54         u_int32_t se,i0,i1;
  55         GET_LDOUBLE_WORDS(se,i0,i1,x);
  56         if(se<0x3fff) {                 /* x < 1 */
  57             return (x-x)/(x-x);
  58         } else if(se >=0x401b) {        /* x > 2**28 */
  59             if(se >=0x7fff) {           /* x is inf of NaN */
  60                 return x+x;
  61             } else
  62                 return __ieee754_logl(x)+ln2;   /* acoshl(huge)=logl(2x) */
  63         } else if(((se-0x3fff)|i0|i1)==0) {
  64             return 0.0;                 /* acosh(1) = 0 */
  65         } else if (se > 0x4000) {       /* 2**28 > x > 2 */
  66             t=x*x;
  67             return __ieee754_logl(2.0*x-one/(x+__ieee754_sqrtl(t-one)));
  68         } else {                        /* 1<x<2 */
  69             t = x-one;
  70             return __log1pl(t+__sqrtl(2.0*t+t*t));
  71         }
  72 }