1/* 2 * Copyright (c) 1998, 2001, Oracle and/or its affiliates. All rights reserved. 3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 4 * 5 * This code is free software; you can redistribute it and/or modify it 6 * under the terms of the GNU General Public License version 2 only, as 7 * published by the Free Software Foundation. Oracle designates this 8 * particular file as subject to the "Classpath" exception as provided 9 * by Oracle in the LICENSE file that accompanied this code. 10 * 11 * This code is distributed in the hope that it will be useful, but WITHOUT 12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14 * version 2 for more details (a copy is included in the LICENSE file that 15 * accompanied this code). 16 * 17 * You should have received a copy of the GNU General Public License version 18 * 2 along with this work; if not, write to the Free Software Foundation, 19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 20 * 21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA 22 * or visit www.oracle.com if you need additional information or have any 23 * questions. 24 */ 25 26/* __ieee754_atanh(x) 27 * Method : 28 * 1.Reduced x to positive by atanh(-x) = -atanh(x) 29 * 2.For x>=0.5 30 * 1 2x x 31 * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) 32 * 2 1 - x 1 - x 33 * 34 * For x<0.5 35 * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) 36 * 37 * Special cases: 38 * atanh(x) is NaN if |x| > 1 with signal; 39 * atanh(NaN) is that NaN with no signal; 40 * atanh(+-1) is +-INF with signal. 41 * 42 */ 43 44#include "fdlibm.h" 45 46#ifdef __STDC__ 47static const double one = 1.0, huge = 1e300; 48#else 49static double one = 1.0, huge = 1e300; 50#endif 51 52static double zero = 0.0; 53 54#ifdef __STDC__ 55 double __ieee754_atanh(double x) 56#else 57 double __ieee754_atanh(x) 58 double x; 59#endif 60{ 61 double t; 62 int hx,ix; 63 unsigned lx; 64 hx = __HI(x); /* high word */ 65 lx = __LO(x); /* low word */ 66 ix = hx&0x7fffffff; 67 if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ 68 return (x-x)/(x-x); 69 if(ix==0x3ff00000) 70 return x/zero; 71 if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ 72 __HI(x) = ix; /* x <- |x| */ 73 if(ix<0x3fe00000) { /* x < 0.5 */ 74 t = x+x; 75 t = 0.5*log1p(t+t*x/(one-x)); 76 } else 77 t = 0.5*log1p((x+x)/(one-x)); 78 if(hx>=0) return t; else return -t; 79} 80