Home | History | Annotate | Download | only in src 
      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 (at) 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 }
     92