Back to index

glibc  2.9
s_asinhl.c
Go to the documentation of this file.
00001 /* @(#)s_asinh.c 5.1 93/09/24 */
00002 /*
00003  * ====================================================
00004  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00005  *
00006  * Developed at SunPro, a Sun Microsystems, Inc. business.
00007  * Permission to use, copy, modify, and distribute this
00008  * software is freely granted, provided that this notice
00009  * is preserved.
00010  * ====================================================
00011  */
00012 
00013 #if defined(LIBM_SCCS) && !defined(lint)
00014 static char rcsid[] = "$NetBSD: s_asinh.c,v 1.9 1995/05/12 04:57:37 jtc Exp $";
00015 #endif
00016 
00017 /* asinh(x)
00018  * Method :
00019  *     Based on
00020  *            asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
00021  *     we have
00022  *     asinh(x) := x  if  1+x*x=1,
00023  *             := sign(x)*(log(x)+ln2)) for large |x|, else
00024  *             := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
00025  *             := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
00026  */
00027 
00028 #include "math.h"
00029 #include "math_private.h"
00030 #include <math_ldbl_opt.h>
00031 
00032 #ifdef __STDC__
00033 static const long double
00034 #else
00035 static long double
00036 #endif
00037 one =  1.00000000000000000000e+00L, /* 0x3ff0000000000000, 0 */
00038 ln2 =  0.6931471805599453094172321214581766L, /* 0x3fe62e42fefa39ef, 0x3c7abc9e3b398040 */
00039 huge=  1.00000000000000000000e+300L;
00040 
00041 #ifdef __STDC__
00042        long double __asinhl(long double x)
00043 #else
00044        long double __asinhl(x)
00045        long double x;
00046 #endif
00047 {
00048        long double t,w;
00049        int64_t hx,ix;
00050        GET_LDOUBLE_MSW64(hx,x);
00051        ix = hx&0x7fffffffffffffffLL;
00052        if(ix>=0x7ff0000000000000LL) return x+x;  /* x is inf or NaN */
00053        if(ix< 0x3e20000000000000LL) {     /* |x|<2**-29 */
00054            if(huge+x>one) return x;       /* return x inexact except 0 */
00055        }
00056        if(ix>0x41b0000000000000LL) {      /* |x| > 2**28 */
00057            w = __ieee754_logl(fabs(x))+ln2;
00058        } else if (ix>0x4000000000000000LL) {     /* 2**28 > |x| > 2.0 */
00059            t = fabs(x);
00060            w = __ieee754_logl(2.0*t+one/(__ieee754_sqrtl(x*x+one)+t));
00061        } else {             /* 2.0 > |x| > 2**-29 */
00062            t = x*x;
00063            w =__log1pl(fabsl(x)+t/(one+__ieee754_sqrtl(one+t)));
00064        }
00065        if(hx>0) return w; else return -w;
00066 }
00067 long_double_symbol (libm, __asinhl, asinhl);