1*1e651e1eSRoland Levillain 2*1e651e1eSRoland Levillain /* @(#)s_asinh.c 1.3 95/01/18 */ 3*1e651e1eSRoland Levillain /* 4*1e651e1eSRoland Levillain * ==================================================== 5*1e651e1eSRoland Levillain * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6*1e651e1eSRoland Levillain * 7*1e651e1eSRoland Levillain * Developed at SunSoft, a Sun Microsystems, Inc. business. 8*1e651e1eSRoland Levillain * Permission to use, copy, modify, and distribute this 9*1e651e1eSRoland Levillain * software is freely granted, provided that this notice 10*1e651e1eSRoland Levillain * is preserved. 11*1e651e1eSRoland Levillain * ==================================================== 12*1e651e1eSRoland Levillain */ 13*1e651e1eSRoland Levillain 14*1e651e1eSRoland Levillain /* asinh(x) 15*1e651e1eSRoland Levillain * Method : 16*1e651e1eSRoland Levillain * Based on 17*1e651e1eSRoland Levillain * asinh(x) = sign(x) * log [ |x| + ieee_sqrt(x*x+1) ] 18*1e651e1eSRoland Levillain * we have 19*1e651e1eSRoland Levillain * asinh(x) := x if 1+x*x=1, 20*1e651e1eSRoland Levillain * := sign(x)*(ieee_log(x)+ln2)) for large |x|, else 21*1e651e1eSRoland Levillain * := sign(x)*ieee_log(2|x|+1/(|x|+ieee_sqrt(x*x+1))) if|x|>2, else 22*1e651e1eSRoland Levillain * := sign(x)*ieee_log1p(|x| + x^2/(1 + ieee_sqrt(1+x^2))) 23*1e651e1eSRoland Levillain */ 24*1e651e1eSRoland Levillain 25*1e651e1eSRoland Levillain #include "fdlibm.h" 26*1e651e1eSRoland Levillain 27*1e651e1eSRoland Levillain #ifdef __STDC__ 28*1e651e1eSRoland Levillain static const double 29*1e651e1eSRoland Levillain #else 30*1e651e1eSRoland Levillain static double 31*1e651e1eSRoland Levillain #endif 32*1e651e1eSRoland Levillain one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 33*1e651e1eSRoland Levillain ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ 34*1e651e1eSRoland Levillain huge= 1.00000000000000000000e+300; 35*1e651e1eSRoland Levillain 36*1e651e1eSRoland Levillain #ifdef __STDC__ asinh(double x)37*1e651e1eSRoland Levillain double asinh(double x) 38*1e651e1eSRoland Levillain #else 39*1e651e1eSRoland Levillain double asinh(x) 40*1e651e1eSRoland Levillain double x; 41*1e651e1eSRoland Levillain #endif 42*1e651e1eSRoland Levillain { 43*1e651e1eSRoland Levillain double t,w; 44*1e651e1eSRoland Levillain int hx,ix; 45*1e651e1eSRoland Levillain hx = __HI(x); 46*1e651e1eSRoland Levillain ix = hx&0x7fffffff; 47*1e651e1eSRoland Levillain if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */ 48*1e651e1eSRoland Levillain if(ix< 0x3e300000) { /* |x|<2**-28 */ 49*1e651e1eSRoland Levillain if(huge+x>one) return x; /* return x inexact except 0 */ 50*1e651e1eSRoland Levillain } 51*1e651e1eSRoland Levillain if(ix>0x41b00000) { /* |x| > 2**28 */ 52*1e651e1eSRoland Levillain w = __ieee754_log(ieee_fabs(x))+ln2; 53*1e651e1eSRoland Levillain } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ 54*1e651e1eSRoland Levillain t = ieee_fabs(x); 55*1e651e1eSRoland Levillain w = __ieee754_log(2.0*t+one/(ieee_sqrt(x*x+one)+t)); 56*1e651e1eSRoland Levillain } else { /* 2.0 > |x| > 2**-28 */ 57*1e651e1eSRoland Levillain t = x*x; 58*1e651e1eSRoland Levillain w =ieee_log1p(ieee_fabs(x)+t/(one+ieee_sqrt(one+t))); 59*1e651e1eSRoland Levillain } 60*1e651e1eSRoland Levillain if(hx>0) return w; else return -w; 61*1e651e1eSRoland Levillain } 62