e_atanh.c revision 141296
1251881Speter 2251881Speter/* @(#)e_atanh.c 1.3 95/01/18 */ 3251881Speter/* 4251881Speter * ==================================================== 5251881Speter * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6251881Speter * 7251881Speter * Developed at SunSoft, a Sun Microsystems, Inc. business. 8251881Speter * Permission to use, copy, modify, and distribute this 9251881Speter * software is freely granted, provided that this notice 10251881Speter * is preserved. 11251881Speter * ==================================================== 12251881Speter * 13251881Speter */ 14251881Speter 15251881Speter#ifndef lint 16251881Speterstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/e_atanh.c 141296 2005-02-04 18:26:06Z das $"; 17251881Speter#endif 18251881Speter 19251881Speter/* __ieee754_atanh(x) 20251881Speter * Method : 21251881Speter * 1.Reduced x to positive by atanh(-x) = -atanh(x) 22251881Speter * 2.For x>=0.5 23251881Speter * 1 2x x 24251881Speter * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) 25251881Speter * 2 1 - x 1 - x 26251881Speter * 27251881Speter * For x<0.5 28251881Speter * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) 29251881Speter * 30251881Speter * Special cases: 31251881Speter * atanh(x) is NaN if |x| > 1 with signal; 32251881Speter * atanh(NaN) is that NaN with no signal; 33251881Speter * atanh(+-1) is +-INF with signal. 34251881Speter * 35251881Speter */ 36251881Speter 37251881Speter#include "math.h" 38251881Speter#include "math_private.h" 39251881Speter 40251881Speterstatic const double one = 1.0, huge = 1e300; 41251881Speterstatic const double zero = 0.0; 42251881Speter 43251881Speterdouble 44251881Speter__ieee754_atanh(double x) 45251881Speter{ 46251881Speter double t; 47251881Speter int32_t hx,ix; 48251881Speter u_int32_t lx; 49251881Speter EXTRACT_WORDS(hx,lx,x); 50251881Speter ix = hx&0x7fffffff; 51251881Speter if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ 52251881Speter return (x-x)/(x-x); 53251881Speter if(ix==0x3ff00000) 54251881Speter return x/zero; 55251881Speter if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ 56251881Speter SET_HIGH_WORD(x,ix); 57251881Speter if(ix<0x3fe00000) { /* x < 0.5 */ 58289180Speter t = x+x; 59251881Speter t = 0.5*log1p(t+t*x/(one-x)); 60251881Speter } else 61251881Speter t = 0.5*log1p((x+x)/(one-x)); 62251881Speter if(hx>=0) return t; else return -t; 63251881Speter} 64251881Speter