POK
/home/jaouen/pok_official/pok/trunk/libpok/libm/e_atanh.c
00001 /*
00002  *                               POK header
00003  * 
00004  * The following file is a part of the POK project. Any modification should
00005  * made according to the POK licence. You CANNOT use this file or a part of
00006  * this file is this part of a file for your own project
00007  *
00008  * For more information on the POK licence, please see our LICENCE FILE
00009  *
00010  * Please follow the coding guidelines described in doc/CODING_GUIDELINES
00011  *
00012  *                                      Copyright (c) 2007-2009 POK team 
00013  *
00014  * Created by julien on Fri Jan 30 14:41:34 2009 
00015  */
00016 
00017 /* @(#)e_atanh.c 5.1 93/09/24 */
00018 /*
00019  * ====================================================
00020  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00021  *
00022  * Developed at SunPro, a Sun Microsystems, Inc. business.
00023  * Permission to use, copy, modify, and distribute this
00024  * software is freely granted, provided that this notice
00025  * is preserved.
00026  * ====================================================
00027  */
00028 
00029 /* __ieee754_atanh(x)
00030  * Method :
00031  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
00032  *    2.For x>=0.5
00033  *                  1              2x                          x
00034  *      atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
00035  *                  2             1 - x                      1 - x
00036  *
00037  *      For x<0.5
00038  *      atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
00039  *
00040  * Special cases:
00041  *      atanh(x) is NaN if |x| > 1 with signal;
00042  *      atanh(NaN) is that NaN with no signal;
00043  *      atanh(+-1) is +-INF with signal.
00044  *
00045  */
00046 #ifdef POK_NEEDS_LIBMATH
00047 
00048 #include <libm.h>
00049 #include "math_private.h"
00050 
00051 static const double one = 1.0, huge = 1e300;
00052 
00053 static const double zero = 0.0;
00054 
00055 double
00056 __ieee754_atanh(double x)
00057 {
00058         double t;
00059         int32_t hx,ix;
00060         uint32_t lx;
00061         EXTRACT_WORDS(hx,lx,x);
00062         ix = hx&0x7fffffff;
00063         if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
00064             return (x-x)/(x-x);
00065         if(ix==0x3ff00000)
00066             return x/zero;
00067         if(ix<0x3e300000&&(huge+x)>zero) return x;      /* x<2**-28 */
00068         SET_HIGH_WORD(x,ix);
00069         if(ix<0x3fe00000) {             /* x < 0.5 */
00070             t = x+x;
00071             t = 0.5*log1p(t+t*x/(one-x));
00072         } else
00073             t = 0.5*log1p((x+x)/(one-x));
00074         if(hx>=0) return t; else return -t;
00075 }
00076 #endif