lib/libm/atanh.c

   1 /*                                                      atanh.c
   2  *
   3  *      Inverse hyperbolic tangent
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, atanh();
  10  *
  11  * y = atanh( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns inverse hyperbolic tangent of argument in the range
  18  * MINLOG to MAXLOG.
  19  *
  20  * If |x| < 0.5, the rational form x + x**3 P(x)/Q(x) is
  21  * employed.  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  *    DEC       -1,1        50000       2.4e-17     6.4e-18
  31  *    IEEE      -1,1        30000       1.9e-16     5.2e-17
  32  *
  33  */
  34 \f
  35 /*                                              atanh.c */
  36
  37
  38 /*
  39 Cephes Math Library Release 2.8:  June, 2000
  40 Copyright (C) 1987, 1995, 2000 by Stephen L. Moshier
  41 */
  42
  43 #include "mconf.h"
  44
  45 #ifdef UNK
  46 static double P[] = {
  47 -8.54074331929669305196E-1,
  48  1.20426861384072379242E1,
  49 -4.61252884198732692637E1,
  50  6.54566728676544377376E1,
  51 -3.09092539379866942570E1
  52 };
  53 static double Q[] = {
  54 /* 1.00000000000000000000E0,*/
  55 -1.95638849376911654834E1,
  56  1.08938092147140262656E2,
  57 -2.49839401325893582852E2,
  58  2.52006675691344555838E2,
  59 -9.27277618139601130017E1
  60 };
  61 #endif
  62 #ifdef DEC
  63 static unsigned short P[] = {
  64 0140132,0122235,0105775,0130300,
  65 0041100,0127327,0124407,0034722,
  66 0141470,0100113,0115607,0130535,
  67 0041602,0164721,0003257,0013673,
  68 0141367,0043046,0166673,0045750
  69 };
  70 static unsigned short Q[] = {
  71 /*0040200,0000000,0000000,0000000,*/
  72 0141234,0101326,0015460,0134564,
  73 0041731,0160115,0116451,0032045,
  74 0142171,0153343,0000532,0167226,
  75 0042174,0000665,0077604,0000310,
  76 0141671,0072235,0031114,0074377
  77 };
  78 #endif
  79
  80 #ifdef IBMPC
  81 static unsigned short P[] = {
  82 0xb618,0xb17f,0x5493,0xbfeb,
  83 0xe73a,0xf520,0x15da,0x4028,
  84 0xf62c,0x7370,0x1009,0xc047,
  85 0xe2f7,0x20d5,0x5d3a,0x4050,
  86 0x697d,0xddb7,0xe8c4,0xc03e
  87 };
  88 static unsigned short Q[] = {
  89 /*0x0000,0x0000,0x0000,0x3ff0,*/
  90 0x172f,0xc366,0x905a,0xc033,
  91 0x2685,0xb3a5,0x3c09,0x405b,
  92 0x5dd3,0x602b,0x3adc,0xc06f,
  93 0x8019,0xaff0,0x8036,0x406f,
  94 0x8f20,0xa649,0x2e93,0xc057
  95 };
  96 #endif
  97
  98 #ifdef MIEEE
  99 static unsigned short P[] = {
 100 0xbfeb,0x5493,0xb17f,0xb618,
 101 0x4028,0x15da,0xf520,0xe73a,
 102 0xc047,0x1009,0x7370,0xf62c,
 103 0x4050,0x5d3a,0x20d5,0xe2f7,
 104 0xc03e,0xe8c4,0xddb7,0x697d
 105 };
 106 static unsigned short Q[] = {
 107 0xc033,0x905a,0xc366,0x172f,
 108 0x405b,0x3c09,0xb3a5,0x2685,
 109 0xc06f,0x3adc,0x602b,0x5dd3,
 110 0x406f,0x8036,0xaff0,0x8019,
 111 0xc057,0x2e93,0xa649,0x8f20
 112 };
 113 #endif
 114
 115 #ifdef ANSIPROT
 116 extern double fabs ( double );
 117 extern double log ( double x );
 118 extern double polevl ( double x, void *P, int N );
 119 extern double p1evl ( double x, void *P, int N );
 120 #else
 121 double fabs(), log(), polevl(), p1evl();
 122 #endif
 123 extern double INFINITY, NAN;
 124
 125 double atanh(x)
 126 double x;
 127 {
 128 double s, z;
 129
 130 #ifdef MINUSZERO
 131 if( x == 0.0 )
 132         return(x);
 133 #endif
 134 z = fabs(x);
 135 if( z >= 1.0 )
 136         {
 137         if( x == 1.0 )
 138                 return( INFINITY );
 139         if( x == -1.0 )
 140                 return( -INFINITY );
 141         mtherr( "atanh", DOMAIN );
 142         return( NAN );
 143         }
 144
 145 if( z < 1.0e-7 )
 146         return(x);
 147
 148 if( z < 0.5 )
 149         {
 150         z = x * x;
 151         s = x   +  x * z * (polevl(z, P, 4) / p1evl(z, Q, 5));
 152         return(s);
 153         }
 154
 155 return( 0.5 * log((1.0+x)/(1.0-x)) );
 156 }