Back to index

glibc  2.9
s_asinh.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 
00031 #ifdef __STDC__
00032 static const double
00033 #else
00034 static double
00035 #endif
00036 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
00037 ln2 =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
00038 huge=  1.00000000000000000000e+300;
00039 
00040 #ifdef __STDC__
00041        double __asinh(double x)
00042 #else
00043        double __asinh(x)
00044        double x;
00045 #endif
00046 {
00047        double t,w;
00048        int32_t hx,ix;
00049        GET_HIGH_WORD(hx,x);
00050        ix = hx&0x7fffffff;
00051        if(ix>=0x7ff00000) return x+x;     /* x is inf or NaN */
00052        if(ix< 0x3e300000) { /* |x|<2**-28 */
00053            if(huge+x>one) return x;       /* return x inexact except 0 */
00054        }
00055        if(ix>0x41b00000) {  /* |x| > 2**28 */
00056            w = __ieee754_log(fabs(x))+ln2;
00057        } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
00058            t = fabs(x);
00059            w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t));
00060        } else {             /* 2.0 > |x| > 2**-28 */
00061            t = x*x;
00062            w =__log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));
00063        }
00064        if(hx>0) return w; else return -w;
00065 }
00066 weak_alias (__asinh, asinh)
00067 #ifdef NO_LONG_DOUBLE
00068 strong_alias (__asinh, __asinhl)
00069 weak_alias (__asinh, asinhl)
00070 #endif