RTAI API: base/math/e_atanh.c Source File

00001 /* @(#)e_atanh.c 5.1 93/09/24 */ 00002 /* 00003 * ==================================================== 00004 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 00005 * 00006 * Developed at SunPro, a Sun Microsystems, Inc. business. 00007 * Permission to use, copy, modify, and distribute this 00008 * software is freely granted, provided that this notice 00009 * is preserved. 00010 * ==================================================== 00011 */ 00012 00013 #if defined(LIBM_SCCS) && !defined(lint) 00014 static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $"; 00015 #endif 00016 00017 /* __ieee754_atanh(x) 00018 * Method : 00019 * 1.Reduced x to positive by atanh(-x) = -atanh(x) 00020 * 2.For x>=0.5 00021 * 1 2x x 00022 * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) 00023 * 2 1 - x 1 - x 00024 * 00025 * For x<0.5 00026 * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) 00027 * 00028 * Special cases: 00029 * atanh(x) is NaN if |x| > 1 with signal; 00030 * atanh(NaN) is that NaN with no signal; 00031 * atanh(+-1) is +-INF with signal. 00032 * 00033 */ 00034 00035 #include "math.h" 00036 #include "mathP.h" 00037 00038 #ifdef __STDC__ 00039 static const double one = 1.0, huge = 1e300; 00040 #else 00041 static double one = 1.0, huge = 1e300; 00042 #endif 00043 00044 #ifdef __STDC__ 00045 static const double zero = 0.0; 00046 #else 00047 static double zero = 0.0; 00048 #endif 00049 00050 #ifdef __STDC__ 00051 double __ieee754_atanh(double x) 00052 #else 00053 double __ieee754_atanh(x) 00054 double x; 00055 #endif 00056 { 00057 double t; 00058 int32_t hx,ix; 00059 u_int32_t lx; 00060 EXTRACT_WORDS(hx,lx,x); 00061 ix = hx&0x7fffffff; 00062 if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ 00063 return (x-x)/(x-x); 00064 if(ix==0x3ff00000) 00065 return x/zero; 00066 if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ 00067 SET_HIGH_WORD(x,ix); 00068 if(ix<0x3fe00000) { /* x < 0.5 */ 00069 t = x+x; 00070 t = 0.5*log1p(t+t*x/(one-x)); 00071 } else 00072 t = 0.5*log1p((x+x)/(one-x)); 00073 if(hx>=0) return t; else return -t; 00074 }