newlib/libm/ld/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
  15 #include <sys/cdefs.h>
  16 __FBSDID("$FreeBSD$");
  17
  18 /*
  19  * See s_asinh.c for complete comments.
  20  *
  21  * Converted to long double by David Schultz <das@FreeBSD.ORG> and
  22  * Bruce D. Evans.
  23  */
  24
  25 #include <float.h>
  26 #ifdef __i386__
  27 #include <ieeefp.h>
  28 #endif
  29
  30 #include "fpmath.h"
  31 #include "math.h"
  32 #include "math_private.h"
  33
  34 /* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */
  35 /* EXP_TINY is the threshold below which we use asinh(x) ~= x. */
  36 #if LDBL_MANT_DIG == 64
  37 #define EXP_LARGE       34
  38 #define EXP_TINY        -34
  39 #elif LDBL_MANT_DIG == 113
  40 #define EXP_LARGE       58
  41 #define EXP_TINY        -58
  42 #else
  43 #error "Unsupported long double format"
  44 #endif
  45
  46 #if LDBL_MAX_EXP != 0x4000
  47 /* We also require the usual expsign encoding. */
  48 #error "Unsupported long double format"
  49 #endif
  50
  51 #define BIAS    (LDBL_MAX_EXP - 1)
  52
  53 static const double
  54 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
  55 huge=  1.00000000000000000000e+300;
  56
  57 #if LDBL_MANT_DIG == 64
  58 static const union IEEEl2bits
  59 u_ln2 =  LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
  60 #define ln2     u_ln2.e
  61 #elif LDBL_MANT_DIG == 113
  62 static const long double
  63 ln2 =  6.93147180559945309417232121458176568e-1L;       /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
  64 #else
  65 #error "Unsupported long double format"
  66 #endif
  67
  68 long double
  69 asinhl(long double x)
  70 {
  71         long double t, w;
  72         uint16_t hx, ix;
  73
  74         ENTERI();
  75         GET_LDBL_EXPSIGN(hx, x);
  76         ix = hx & 0x7fff;
  77         if (ix >= 0x7fff) RETURNI(x+x); /* x is inf, NaN or misnormal */
  78         if (ix < BIAS + EXP_TINY) {     /* |x| < TINY, or misnormal */
  79             if (huge + x > one) RETURNI(x);     /* return x inexact except 0 */
  80         }
  81         if (ix >= BIAS + EXP_LARGE) {   /* |x| >= LARGE, or misnormal */
  82             w = logl(fabsl(x))+ln2;
  83         } else if (ix >= 0x4000) {      /* LARGE > |x| >= 2.0, or misnormal */
  84             t = fabsl(x);
  85             w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
  86         } else {                /* 2.0 > |x| >= TINY, or misnormal */
  87             t = x*x;
  88             w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
  89         }
  90         RETURNI((hx & 0x8000) == 0 ? w : -w);
  91 }