lib/libm/src/s_asinhl.c

   1 /* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
   2
   3 /* @(#)s_asinh.c 5.1 93/09/24 */
   4 /*
   5  * ====================================================
   6  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   7  *
   8  * Developed at SunPro, a Sun Microsystems, Inc. business.
   9  * Permission to use, copy, modify, and distribute this
  10  * software is freely granted, provided that this notice
  11  * is preserved.
  12  * ====================================================
  13  *
  14  * $FreeBSD: head/lib/msun/src/s_asinhl.c 251599 2013-06-10 06:04:58Z das $
  15  */
  16
  17 /*
  18  * See s_asinh.c for complete comments.
  19  *
  20  * Converted to long double by David Schultz <das@FreeBSD.ORG> and
  21  * Bruce D. Evans.
  22  */
  23
  24 #include <float.h>
  25 #ifdef __i386__
  26 #include <ieeefp.h>
  27 #endif
  28
  29 #include "fpmath.h"
  30 #include "math.h"
  31 #include "math_private.h"
  32
  33 /* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */
  34 /* EXP_TINY is the threshold below which we use asinh(x) ~= x. */
  35 #if LDBL_MANT_DIG == 64
  36 #define EXP_LARGE       34
  37 #define EXP_TINY        -34
  38 #elif LDBL_MANT_DIG == 113
  39 #define EXP_LARGE       58
  40 #define EXP_TINY        -58
  41 #else
  42 #error "Unsupported long double format"
  43 #endif
  44
  45 #if LDBL_MAX_EXP != 0x4000
  46 /* We also require the usual expsign encoding. */
  47 #error "Unsupported long double format"
  48 #endif
  49
  50 #define BIAS    (LDBL_MAX_EXP - 1)
  51
  52 static const double
  53 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
  54 huge=  1.00000000000000000000e+300;
  55
  56 #if LDBL_MANT_DIG == 64
  57 static const union IEEEl2bits
  58 u_ln2 =  LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
  59 #define ln2     u_ln2.e
  60 #elif LDBL_MANT_DIG == 113
  61 static const long double
  62 ln2 =  6.93147180559945309417232121458176568e-1L;       /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
  63 #else
  64 #error "Unsupported long double format"
  65 #endif
  66
  67 long double
  68 asinhl(long double x)
  69 {
  70         long double t, w;
  71         uint16_t hx, ix;
  72
  73         ENTERI();
  74         GET_LDBL_EXPSIGN(hx, x);
  75         ix = hx & 0x7fff;
  76         if (ix >= 0x7fff) RETURNI(x+x); /* x is inf, NaN or misnormal */
  77         if (ix < BIAS + EXP_TINY) {     /* |x| < TINY, or misnormal */
  78             if (huge + x > one) RETURNI(x);     /* return x inexact except 0 */
  79         }
  80         if (ix >= BIAS + EXP_LARGE) {   /* |x| >= LARGE, or misnormal */
  81             w = logl(fabsl(x))+ln2;
  82         } else if (ix >= 0x4000) {      /* LARGE > |x| >= 2.0, or misnormal */
  83             t = fabsl(x);
  84             w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
  85         } else {                /* 2.0 > |x| >= TINY, or misnormal */
  86             t = x*x;
  87             w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
  88         }
  89         RETURNI((hx & 0x8000) == 0 ? w : -w);
  90 }