Back to index

glibc  2.9
e_atanhl.c
Go to the documentation of this file.
00001 /* s_atanhl.c -- long double version of s_atan.c.
00002  * Conversion to long double by Ulrich Drepper,
00003  * Cygnus Support, drepper@cygnus.com.
00004  */
00005 
00006 /*
00007  * ====================================================
00008  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00009  *
00010  * Developed at SunPro, a Sun Microsystems, Inc. business.
00011  * Permission to use, copy, modify, and distribute this
00012  * software is freely granted, provided that this notice
00013  * is preserved.
00014  * ====================================================
00015  */
00016 
00017 #if defined(LIBM_SCCS) && !defined(lint)
00018 static char rcsid[] = "$NetBSD: $";
00019 #endif
00020 
00021 /* __ieee754_atanhl(x)
00022  * Method :
00023  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
00024  *    2.For x>=0.5
00025  *                   1              2x                          x
00026  *     atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
00027  *                   2             1 - x                      1 - x
00028  *
00029  *     For x<0.5
00030  *     atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
00031  *
00032  * Special cases:
00033  *     atanhl(x) is NaN if |x| > 1 with signal;
00034  *     atanhl(NaN) is that NaN with no signal;
00035  *     atanhl(+-1) is +-INF with signal.
00036  *
00037  */
00038 
00039 #include "math.h"
00040 #include "math_private.h"
00041 
00042 #ifdef __STDC__
00043 static const long double one = 1.0, huge = 1e4900L;
00044 #else
00045 static long double one = 1.0, huge = 1e4900L;
00046 #endif
00047 
00048 #ifdef __STDC__
00049 static const long double zero = 0.0;
00050 #else
00051 static double long zero = 0.0;
00052 #endif
00053 
00054 #ifdef __STDC__
00055        long double __ieee754_atanhl(long double x)
00056 #else
00057        long double __ieee754_atanhl(x)
00058        long double x;
00059 #endif
00060 {
00061        long double t;
00062        int32_t ix;
00063        u_int32_t se,i0,i1;
00064        GET_LDOUBLE_WORDS(se,i0,i1,x);
00065        ix = se&0x7fff;
00066        if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31))>0x3fff)
00067          /* |x|>1 */
00068            return (x-x)/(x-x);
00069        if(ix==0x3fff)
00070            return x/zero;
00071        if(ix<0x3fe3&&(huge+x)>zero) return x;    /* x<2**-28 */
00072        SET_LDOUBLE_EXP(x,ix);
00073        if(ix<0x3ffe) {             /* x < 0.5 */
00074            t = x+x;
00075            t = 0.5*__log1pl(t+t*x/(one-x));
00076        } else
00077            t = 0.5*__log1pl((x+x)/(one-x));
00078        if(se<=0x7fff) return t; else return -t;
00079 }