newlib/libm/mathfp/e_atanh.c

   1
   2 /* @(#)e_atanh.c 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  *
  13  */
  14
  15 /*
  16 FUNCTION
  17         <<atanh>>, <<atanhf>>---inverse hyperbolic tangent
  18
  19 INDEX
  20         atanh
  21 INDEX
  22         atanhf
  23
  24 SYNOPSIS
  25         #include <math.h>
  26         double atanh(double <[x]>);
  27         float atanhf(float <[x]>);
  28
  29 DESCRIPTION
  30         <<atanh>> calculates the inverse hyperbolic tangent of <[x]>.
  31
  32         <<atanhf>> is identical, other than taking and returning
  33         <<float>> values.
  34
  35 RETURNS
  36         <<atanh>> and <<atanhf>> return the calculated value.
  37
  38         If
  39         @ifnottex
  40         |<[x]>|
  41         @end ifnottex
  42         @tex
  43         $|x|$
  44         @end tex
  45         is greater than 1, the global <<errno>> is set to <<EDOM>> and
  46         the result is a NaN.  A <<DOMAIN error>> is reported.
  47
  48         If
  49         @ifnottex
  50         |<[x]>|
  51         @end ifnottex
  52         @tex
  53         $|x|$
  54         @end tex
  55         is 1, the global <<errno>> is set to <<EDOM>>; and the result is
  56         infinity with the same sign as <<x>>.  A <<SING error>> is reported.
  57
  58 PORTABILITY
  59         Neither <<atanh>> nor <<atanhf>> are ANSI C.
  60
  61 QUICKREF
  62         atanh - pure
  63         atanhf - pure
  64
  65
  66 */
  67
  68 /* atanh(x)
  69  * Method :
  70  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
  71  *    2.For x>=0.5
  72  *                  1              2x                          x
  73  *      atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
  74  *                  2             1 - x                      1 - x
  75  *
  76  *      For x<0.5
  77  *      atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
  78  *
  79  * Special cases:
  80  *      atanh(x) is NaN if |x| > 1 with signal;
  81  *      atanh(NaN) is that NaN with no signal;
  82  *      atanh(+-1) is +-INF with signal.
  83  *
  84  */
  85
  86 #include "fdlibm.h"
  87
  88 #ifndef _DOUBLE_IS_32BITS
  89
  90 #ifdef __STDC__
  91 static const double one = 1.0, huge = 1e300;
  92 #else
  93 static double one = 1.0, huge = 1e300;
  94 #endif
  95
  96 #ifdef __STDC__
  97 static const double zero = 0.0;
  98 #else
  99 static double zero = 0.0;
 100 #endif
 101
 102 #ifdef __STDC__
 103         double atanh(double x)
 104 #else
 105         double atanh(x)
 106         double x;
 107 #endif
 108 {
 109         double t;
 110         __int32_t hx,ix;
 111         __uint32_t lx;
 112         EXTRACT_WORDS(hx,lx,x);
 113         ix = hx&0x7fffffff;
 114         if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
 115             return (x-x)/(x-x);
 116         if(ix==0x3ff00000)
 117             return x/zero;
 118         if(ix<0x3e300000&&(huge+x)>zero) return x;      /* x<2**-28 */
 119         SET_HIGH_WORD(x,ix);
 120         if(ix<0x3fe00000) {             /* x < 0.5 */
 121             t = x+x;
 122             t = 0.5*log1p(t+t*x/(one-x));
 123         } else
 124             t = 0.5*log1p((x+x)/(one-x));
 125         if(hx>=0) return t; else return -t;
 126 }
 127
 128 #endif /* defined(_DOUBLE_IS_32BITS) */