e_atanhl.c revision 22993
1264790Sbapt/* @(#)e_atanh.c 5.1 93/09/24 */ 2264790Sbapt/* 3264790Sbapt * ==================================================== 4264790Sbapt * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5264790Sbapt * 6264790Sbapt * Developed at SunPro, a Sun Microsystems, Inc. business. 7264790Sbapt * Permission to use, copy, modify, and distribute this 8264790Sbapt * software is freely granted, provided that this notice 9264790Sbapt * is preserved. 10264790Sbapt * ==================================================== 11264790Sbapt */ 12264790Sbapt 13264790Sbapt#ifndef lint 14264790Sbaptstatic char rcsid[] = "$Id$"; 15264790Sbapt#endif 16264790Sbapt 17264790Sbapt/* __ieee754_atanh(x) 18264790Sbapt * Method : 19264790Sbapt * 1.Reduced x to positive by atanh(-x) = -atanh(x) 20264790Sbapt * 2.For x>=0.5 21264790Sbapt * 1 2x x 22264790Sbapt * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) 23264790Sbapt * 2 1 - x 1 - x 24264790Sbapt * 25264790Sbapt * For x<0.5 26264790Sbapt * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) 27264790Sbapt * 28264790Sbapt * Special cases: 29264790Sbapt * atanh(x) is NaN if |x| > 1 with signal; 30264790Sbapt * atanh(NaN) is that NaN with no signal; 31264790Sbapt * atanh(+-1) is +-INF with signal. 32264790Sbapt * 33264790Sbapt */ 34264790Sbapt 35264790Sbapt#include "math.h" 36264790Sbapt#include "math_private.h" 37264790Sbapt 38264790Sbapt#ifdef __STDC__ 39264790Sbaptstatic const double one = 1.0, huge = 1e300; 40264790Sbapt#else 41264790Sbaptstatic double one = 1.0, huge = 1e300; 42264790Sbapt#endif 43264790Sbapt 44264790Sbapt#ifdef __STDC__ 45264790Sbaptstatic const double zero = 0.0; 46264790Sbapt#else 47264790Sbaptstatic double zero = 0.0; 48264790Sbapt#endif 49264790Sbapt 50264790Sbapt#ifdef __STDC__ 51264790Sbapt double __ieee754_atanh(double x) 52264790Sbapt#else 53264790Sbapt double __ieee754_atanh(x) 54264790Sbapt double x; 55264790Sbapt#endif 56264790Sbapt{ 57264790Sbapt double t; 58264790Sbapt int32_t hx,ix; 59264790Sbapt u_int32_t lx; 60264790Sbapt EXTRACT_WORDS(hx,lx,x); 61264790Sbapt ix = hx&0x7fffffff; 62264790Sbapt if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ 63264790Sbapt return (x-x)/(x-x); 64264790Sbapt if(ix==0x3ff00000) 65264790Sbapt return x/zero; 66264790Sbapt if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ 67264790Sbapt SET_HIGH_WORD(x,ix); 68264790Sbapt if(ix<0x3fe00000) { /* x < 0.5 */ 69264790Sbapt t = x+x; 70264790Sbapt t = 0.5*log1p(t+t*x/(one-x)); 71264790Sbapt } else 72264790Sbapt t = 0.5*log1p((x+x)/(one-x)); 73264790Sbapt if(hx>=0) return t; else return -t; 74264790Sbapt} 75264790Sbapt