lib/libm/atanhf.c

   1 /*                                                      atanhf.c
   2  *
   3  *      Inverse hyperbolic tangent
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, atanhf();
  10  *
  11  * y = atanhf( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns inverse hyperbolic tangent of argument in the range
  18  * MINLOGF to MAXLOGF.
  19  *
  20  * If |x| < 0.5, a polynomial approximation is used.
  21  * Otherwise,
  22  *        atanh(x) = 0.5 * log( (1+x)/(1-x) ).
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *                      Relative error:
  29  * arithmetic   domain     # trials      peak         rms
  30  *    IEEE      -1,1        100000      1.4e-7      3.1e-8
  31  *
  32  */
  33 \f
  34 /*                                              atanh.c */
  35
  36
  37 /*
  38 Cephes Math Library Release 2.2:  June, 1992
  39 Copyright (C) 1987, 1992 by Stephen L. Moshier
  40 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  41 */
  42
  43 /* Single precision inverse hyperbolic tangent
  44  * test interval: [-0.5, +0.5]
  45  * trials: 10000
  46  * peak relative error: 8.2e-8
  47  * rms relative error: 3.0e-8
  48  */
  49 #include "mconf.h"
  50 extern float MAXNUMF;
  51
  52 #ifdef ANSIC
  53 float logf( float );
  54
  55 float atanhf( float xx )
  56 #else
  57 float logf();
  58
  59 float atanhf(xx)
  60 double xx;
  61 #endif
  62 {
  63 float x, z;
  64
  65 x = xx;
  66 if( x < 0 )
  67         z = -x;
  68 else
  69         z = x;
  70 if( z >= 1.0 )
  71         {
  72         if( x == 1.0 )
  73                 return( MAXNUMF );
  74         if( x == -1.0 )
  75                 return( -MAXNUMF );
  76         mtherr( "atanhl", DOMAIN );
  77         return( MAXNUMF );
  78         }
  79
  80 if( z < 1.0e-4 )
  81         return(x);
  82
  83 if( z < 0.5 )
  84         {
  85         z = x * x;
  86         z =
  87         (((( 1.81740078349E-1 * z
  88           + 8.24370301058E-2) * z
  89           + 1.46691431730E-1) * z
  90           + 1.99782164500E-1) * z
  91           + 3.33337300303E-1) * z * x
  92           + x;
  93         }
  94 else
  95         {
  96         z = 0.5 * logf( (1.0+x)/(1.0-x) );
  97         }
  98 return( z );
  99 }