nsz Git - libm/blob - src/math/acoshl.c

   1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_acoshl.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /* acoshl(x)
  13  * Method :
  14  *      Based on
  15  *              acoshl(x) = logl [ x + sqrtl(x*x-1) ]
  16  *      we have
  17  *              acoshl(x) := logl(x)+ln2,       if x is large; else
  18  *              acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
  19  *              acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
  20  *
  21  * Special cases:
  22  *      acoshl(x) is NaN with signal if x<1.
  23  *      acoshl(NaN) is NaN without signal.
  24  */
  25
  26 #include "libm.h"
  27
  28 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
  29 long double acoshl(long double x)
  30 {
  31         return acosh(x);
  32 }
  33 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
  34 static const long double
  35 one = 1.0,
  36 ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
  37
  38 long double acoshl(long double x)
  39 {
  40         long double t;
  41         uint32_t se,i0,i1;
  42
  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 or NaN */
  48                         return x+x;
  49                 return logl(x) + ln2;  /* acoshl(huge) = logl(2x) */
  50         } else if (((se-0x3fff)|i0|i1) == 0) {
  51                 return 0.0;            /* acosh(1) = 0 */
  52         } else if (se > 0x4000) {  /* x > 2 */
  53                 t = x*x;
  54                 return logl(2.0*x - one/(x + sqrtl(t - one)));
  55         }
  56         /* 1 < x <= 2 */
  57         t = x - one;
  58         return log1pl(t + sqrtl(2.0*t + t*t));
  59 }
  60 #endif