s_asinh.c revision 302408
1185380Ssam/* @(#)s_asinh.c 5.1 93/09/24 */ 2185380Ssam/* 3185380Ssam * ==================================================== 4185380Ssam * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5185380Ssam * 6185380Ssam * Developed at SunPro, a Sun Microsystems, Inc. business. 7185380Ssam * Permission to use, copy, modify, and distribute this 8185380Ssam * software is freely granted, provided that this notice 9185380Ssam * is preserved. 10185380Ssam * ==================================================== 11185380Ssam */ 12185380Ssam 13185380Ssam#include <sys/cdefs.h> 14185380Ssam__FBSDID("$FreeBSD: stable/11/lib/msun/src/s_asinh.c 251599 2013-06-10 06:04:58Z das $"); 15185380Ssam 16185380Ssam/* asinh(x) 17203158Srpaulo * Method : 18185380Ssam * Based on 19185380Ssam * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] 20185380Ssam * we have 21185380Ssam * asinh(x) := x if 1+x*x=1, 22185380Ssam * := sign(x)*(log(x)+ln2)) for large |x|, else 23185380Ssam * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else 24185380Ssam * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) 25185380Ssam */ 26185380Ssam 27185380Ssam#include <float.h> 28185380Ssam 29185380Ssam#include "math.h" 30185380Ssam#include "math_private.h" 31185380Ssam 32185380Ssamstatic const double 33185380Ssamone = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 34185380Ssamln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ 35185380Ssamhuge= 1.00000000000000000000e+300; 36185380Ssam 37185380Ssamdouble 38185380Ssamasinh(double x) 39185380Ssam{ 40185380Ssam double t,w; 41185380Ssam int32_t hx,ix; 42185380Ssam GET_HIGH_WORD(hx,x); 43185380Ssam ix = hx&0x7fffffff; 44185380Ssam if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */ 45185380Ssam if(ix< 0x3e300000) { /* |x|<2**-28 */ 46185380Ssam if(huge+x>one) return x; /* return x inexact except 0 */ 47185380Ssam } 48185380Ssam if(ix>0x41b00000) { /* |x| > 2**28 */ 49185380Ssam w = __ieee754_log(fabs(x))+ln2; 50185380Ssam } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ 51185380Ssam t = fabs(x); 52185380Ssam w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t)); 53185380Ssam } else { /* 2.0 > |x| > 2**-28 */ 54185380Ssam t = x*x; 55185380Ssam w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t))); 56185380Ssam } 57185380Ssam if(hx>0) return w; else return -w; 58185380Ssam} 59185380Ssam 60185380Ssam#if LDBL_MANT_DIG == 53 61185380Ssam__weak_reference(asinh, asinhl); 62185380Ssam#endif 63185380Ssam