Home | History | Annotate | Line # | Download | only in src 
      1 /* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */
      2 
      3 /*
      4  * ====================================================
      5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      6  *
      7  * Developed at SunSoft, a Sun Microsystems, Inc. business.
      8  * Permission to use, copy, modify, and distribute this
      9  * software is freely granted, provided that this notice
     10  * is preserved.
     11  * ====================================================
     12  *
     13  */
     14 
     15 #include <sys/cdefs.h>
     16 
     17 #include "namespace.h"
     18 
     19 #include <float.h>
     20 #include <machine/ieee.h>
     21 
     22 #include "math.h"
     23 #include "math_private.h"
     24 
     25 __weak_alias(atanhl, _atanhl)
     26 
     27 #ifdef __HAVE_LONG_DOUBLE
     28 /*
     29  * See e_atanh.c for complete comments.
     30  *
     31  * Converted to long double by David Schultz <das (at) FreeBSD.ORG> and
     32  * Bruce D. Evans.
     33  */
     34 
     35 /* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
     36 #if LDBL_MANT_DIG == 64
     37 #define	EXP_TINY	-34
     38 #elif LDBL_MANT_DIG == 113
     39 #define	EXP_TINY	-58
     40 #else
     41 #error "Unsupported long double format"
     42 #endif
     43 
     44 #if LDBL_MAX_EXP != 0x4000
     45 /* We also require the usual expsign encoding. */
     46 #error "Unsupported long double format"
     47 #endif
     48 
     49 #define	BIAS	(LDBL_MAX_EXP - 1)
     50 
     51 static const double one = 1.0, huge = 1e300;
     52 static const double zero = 0.0;
     53 
     54 long double
     55 atanhl(long double x)
     56 {
     57 	long double t;
     58 	uint16_t hx, ix;
     59 
     60 	ENTERI();
     61 	GET_LDBL_EXPSIGN(hx, x);
     62 	ix = hx & 0x7fff;
     63 	if (ix >= 0x3fff)		/* |x| >= 1, or NaN or misnormal */
     64 	    RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
     65 	if (ix < BIAS + EXP_TINY && (huge + x) > zero)
     66 	    RETURNI(x);			/* x is tiny */
     67 	SET_LDBL_EXPSIGN(x, ix);
     68 	if (ix < 0x3ffe) {		/* |x| < 0.5, or misnormal */
     69 	    t = x+x;
     70 	    t = 0.5*log1pl(t+t*x/(one-x));
     71 	} else
     72 	    t = 0.5*log1pl((x+x)/(one-x));
     73 	RETURNI((hx & 0x8000) == 0 ? t : -t);
     74 }
     75 #else
     76 long double
     77 atanhl(long double x)
     78 {
     79 	return atanh(x);
     80 }
     81 #endif
     82