Annotation of src/lib/libm/src/s_asinh.c, Revision 1.11
1.1 jtc 1: /* @(#)s_asinh.c 5.1 93/09/24 */
2: /*
3: * ====================================================
4: * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5: *
6: * Developed at SunPro, a Sun Microsystems, Inc. business.
7: * Permission to use, copy, modify, and distribute this
1.11 ! simonb 8: * software is freely granted, provided that this notice
1.1 jtc 9: * is preserved.
10: * ====================================================
11: */
1.3 jtc 12:
1.10 lukem 13: #include <sys/cdefs.h>
1.7 jtc 14: #if defined(LIBM_SCCS) && !defined(lint)
1.11 ! simonb 15: __RCSID("$NetBSD: s_asinh.c,v 1.10 1997/10/09 11:30:50 lukem Exp $");
1.3 jtc 16: #endif
1.1 jtc 17:
18: /* asinh(x)
19: * Method :
1.11 ! simonb 20: * Based on
1.1 jtc 21: * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
22: * we have
23: * asinh(x) := x if 1+x*x=1,
24: * := sign(x)*(log(x)+ln2)) for large |x|, else
25: * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
1.11 ! simonb 26: * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
1.1 jtc 27: */
28:
1.5 jtc 29: #include "math.h"
30: #include "math_private.h"
1.1 jtc 31:
32: #ifdef __STDC__
1.11 ! simonb 33: static const double
1.1 jtc 34: #else
1.11 ! simonb 35: static double
1.1 jtc 36: #endif
37: one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
38: ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
1.11 ! simonb 39: huge= 1.00000000000000000000e+300;
1.1 jtc 40:
41: #ifdef __STDC__
42: double asinh(double x)
43: #else
44: double asinh(x)
45: double x;
46: #endif
1.11 ! simonb 47: {
1.1 jtc 48: double t,w;
1.6 jtc 49: int32_t hx,ix;
1.5 jtc 50: GET_HIGH_WORD(hx,x);
1.1 jtc 51: ix = hx&0x7fffffff;
52: if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */
53: if(ix< 0x3e300000) { /* |x|<2**-28 */
54: if(huge+x>one) return x; /* return x inexact except 0 */
1.11 ! simonb 55: }
1.1 jtc 56: if(ix>0x41b00000) { /* |x| > 2**28 */
57: w = __ieee754_log(fabs(x))+ln2;
58: } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
59: t = fabs(x);
1.9 jtc 60: w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t));
1.1 jtc 61: } else { /* 2.0 > |x| > 2**-28 */
62: t = x*x;
1.9 jtc 63: w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));
1.1 jtc 64: }
65: if(hx>0) return w; else return -w;
66: }
CVSweb <webmaster@jp.NetBSD.org>