1141296Sdas 2141296Sdas/* @(#)e_atanh.c 1.3 95/01/18 */ 32116Sjkh/* 42116Sjkh * ==================================================== 52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 62116Sjkh * 7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 82116Sjkh * Permission to use, copy, modify, and distribute this 9141296Sdas * software is freely granted, provided that this notice 102116Sjkh * is preserved. 112116Sjkh * ==================================================== 12141296Sdas * 132116Sjkh */ 142116Sjkh 15176451Sdas#include <sys/cdefs.h> 16176451Sdas__FBSDID("$FreeBSD$"); 172116Sjkh 182116Sjkh/* __ieee754_atanh(x) 192116Sjkh * Method : 202116Sjkh * 1.Reduced x to positive by atanh(-x) = -atanh(x) 212116Sjkh * 2.For x>=0.5 222116Sjkh * 1 2x x 232116Sjkh * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) 242116Sjkh * 2 1 - x 1 - x 25141296Sdas * 262116Sjkh * For x<0.5 272116Sjkh * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) 282116Sjkh * 292116Sjkh * Special cases: 302116Sjkh * atanh(x) is NaN if |x| > 1 with signal; 312116Sjkh * atanh(NaN) is that NaN with no signal; 322116Sjkh * atanh(+-1) is +-INF with signal. 332116Sjkh * 342116Sjkh */ 352116Sjkh 36251599Sdas#include <float.h> 37251599Sdas 382116Sjkh#include "math.h" 392116Sjkh#include "math_private.h" 402116Sjkh 412116Sjkhstatic const double one = 1.0, huge = 1e300; 422116Sjkhstatic const double zero = 0.0; 432116Sjkh 4497407Salfreddouble 4597407Salfred__ieee754_atanh(double x) 462116Sjkh{ 472116Sjkh double t; 482116Sjkh int32_t hx,ix; 492116Sjkh u_int32_t lx; 502116Sjkh EXTRACT_WORDS(hx,lx,x); 512116Sjkh ix = hx&0x7fffffff; 522116Sjkh if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ 532116Sjkh return (x-x)/(x-x); 54141296Sdas if(ix==0x3ff00000) 552116Sjkh return x/zero; 562116Sjkh if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ 572116Sjkh SET_HIGH_WORD(x,ix); 582116Sjkh if(ix<0x3fe00000) { /* x < 0.5 */ 592116Sjkh t = x+x; 602116Sjkh t = 0.5*log1p(t+t*x/(one-x)); 61141296Sdas } else 622116Sjkh t = 0.5*log1p((x+x)/(one-x)); 632116Sjkh if(hx>=0) return t; else return -t; 642116Sjkh} 65251599Sdas 66251599Sdas#if LDBL_MANT_DIG == 53 67251599Sdas__weak_reference(atanh, atanhl); 68251599Sdas#endif 69