atanh.c 3.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127
  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. /* atanh.c */
  35. /*
  36. Cephes Math Library Release 2.8: June, 2000
  37. Copyright (C) 1987, 1995, 2000 by Stephen L. Moshier
  38. */
  39. #include "mconf.h"
  40. #ifdef UNK
  41. static double P[] = {-8.54074331929669305196E-1, 1.20426861384072379242E1,
  42. -4.61252884198732692637E1, 6.54566728676544377376E1,
  43. -3.09092539379866942570E1};
  44. static double Q[] = {
  45. /* 1.00000000000000000000E0,*/
  46. -1.95638849376911654834E1, 1.08938092147140262656E2,
  47. -2.49839401325893582852E2, 2.52006675691344555838E2,
  48. -9.27277618139601130017E1};
  49. #endif
  50. #ifdef DEC
  51. static unsigned short P[] = {0140132, 0122235, 0105775, 0130300, 0041100,
  52. 0127327, 0124407, 0034722, 0141470, 0100113,
  53. 0115607, 0130535, 0041602, 0164721, 0003257,
  54. 0013673, 0141367, 0043046, 0166673, 0045750};
  55. static unsigned short Q[] = {
  56. /*0040200,0000000,0000000,0000000,*/
  57. 0141234, 0101326, 0015460, 0134564, 0041731, 0160115, 0116451,
  58. 0032045, 0142171, 0153343, 0000532, 0167226, 0042174, 0000665,
  59. 0077604, 0000310, 0141671, 0072235, 0031114, 0074377};
  60. #endif
  61. #ifdef IBMPC
  62. static unsigned short P[] = {0xb618, 0xb17f, 0x5493, 0xbfeb, 0xe73a,
  63. 0xf520, 0x15da, 0x4028, 0xf62c, 0x7370,
  64. 0x1009, 0xc047, 0xe2f7, 0x20d5, 0x5d3a,
  65. 0x4050, 0x697d, 0xddb7, 0xe8c4, 0xc03e};
  66. static unsigned short Q[] = {
  67. /*0x0000,0x0000,0x0000,0x3ff0,*/
  68. 0x172f, 0xc366, 0x905a, 0xc033, 0x2685, 0xb3a5, 0x3c09,
  69. 0x405b, 0x5dd3, 0x602b, 0x3adc, 0xc06f, 0x8019, 0xaff0,
  70. 0x8036, 0x406f, 0x8f20, 0xa649, 0x2e93, 0xc057};
  71. #endif
  72. #ifdef MIEEE
  73. static unsigned short P[] = {0xbfeb, 0x5493, 0xb17f, 0xb618, 0x4028,
  74. 0x15da, 0xf520, 0xe73a, 0xc047, 0x1009,
  75. 0x7370, 0xf62c, 0x4050, 0x5d3a, 0x20d5,
  76. 0xe2f7, 0xc03e, 0xe8c4, 0xddb7, 0x697d};
  77. static unsigned short Q[] = {0xc033, 0x905a, 0xc366, 0x172f, 0x405b,
  78. 0x3c09, 0xb3a5, 0x2685, 0xc06f, 0x3adc,
  79. 0x602b, 0x5dd3, 0x406f, 0x8036, 0xaff0,
  80. 0x8019, 0xc057, 0x2e93, 0xa649, 0x8f20};
  81. #endif
  82. #ifdef ANSIPROT
  83. extern double fabs(double);
  84. extern double log(double x);
  85. extern double polevl(double x, void *P, int N);
  86. extern double p1evl(double x, void *P, int N);
  87. #else
  88. double fabs(), log(), polevl(), p1evl();
  89. #endif
  90. extern double INFINITY, NAN;
  91. double atanh(x) double x;
  92. {
  93. double s, z;
  94. #ifdef MINUSZERO
  95. if (x == 0.0)
  96. return (x);
  97. #endif
  98. z = fabs(x);
  99. if (z >= 1.0) {
  100. if (x == 1.0)
  101. return (INFINITY);
  102. if (x == -1.0)
  103. return (-INFINITY);
  104. mtherr("atanh", DOMAIN);
  105. return (NAN);
  106. }
  107. if (z < 1.0e-7)
  108. return (x);
  109. if (z < 0.5) {
  110. z = x * x;
  111. s = x + x * z * (polevl(z, P, 4) / p1evl(z, Q, 5));
  112. return (s);
  113. }
  114. return (0.5 * log((1.0 + x) / (1.0 - x)));
  115. }