tanh.c 2.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141
  1. /* tanh.c
  2. *
  3. * Hyperbolic tangent
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, tanh();
  10. *
  11. * y = tanh( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns hyperbolic tangent of argument in the range MINLOG to
  18. * MAXLOG.
  19. *
  20. * A rational function is used for |x| < 0.625. The form
  21. * x + x**3 P(x)/Q(x) of Cody _& Waite is employed.
  22. * Otherwise,
  23. * tanh(x) = sinh(x)/cosh(x) = 1 - 2/(exp(2x) + 1).
  24. *
  25. *
  26. *
  27. * ACCURACY:
  28. *
  29. * Relative error:
  30. * arithmetic domain # trials peak rms
  31. * DEC -2,2 50000 3.3e-17 6.4e-18
  32. * IEEE -2,2 30000 2.5e-16 5.8e-17
  33. *
  34. */
  35. /*
  36. Cephes Math Library Release 2.8: June, 2000
  37. Copyright 1984, 1995, 2000 by Stephen L. Moshier
  38. */
  39. #include <math.h>
  40. #ifdef UNK
  41. static double P[] = {
  42. -9.64399179425052238628E-1,
  43. -9.92877231001918586564E1,
  44. -1.61468768441708447952E3
  45. };
  46. static double Q[] = {
  47. /* 1.00000000000000000000E0,*/
  48. 1.12811678491632931402E2,
  49. 2.23548839060100448583E3,
  50. 4.84406305325125486048E3
  51. };
  52. #endif
  53. #ifdef DEC
  54. static unsigned short P[] = {
  55. 0140166,0161335,0053753,0075126,
  56. 0141706,0111520,0070463,0040552,
  57. 0142711,0153001,0101300,0025430
  58. };
  59. static unsigned short Q[] = {
  60. /*0040200,0000000,0000000,0000000,*/
  61. 0041741,0117624,0051300,0156060,
  62. 0043013,0133720,0071251,0127717,
  63. 0043227,0060201,0021020,0020136
  64. };
  65. #endif
  66. #ifdef IBMPC
  67. static unsigned short P[] = {
  68. 0x6f4b,0xaafd,0xdc5b,0xbfee,
  69. 0x682d,0x0e26,0xd26a,0xc058,
  70. 0x0563,0x3058,0x3ac0,0xc099
  71. };
  72. static unsigned short Q[] = {
  73. /*0x0000,0x0000,0x0000,0x3ff0,*/
  74. 0x1b86,0x8a58,0x33f2,0x405c,
  75. 0x35fa,0x0e55,0x76fa,0x40a1,
  76. 0x040c,0x2442,0xec10,0x40b2
  77. };
  78. #endif
  79. #ifdef MIEEE
  80. static unsigned short P[] = {
  81. 0xbfee,0xdc5b,0xaafd,0x6f4b,
  82. 0xc058,0xd26a,0x0e26,0x682d,
  83. 0xc099,0x3ac0,0x3058,0x0563
  84. };
  85. static unsigned short Q[] = {
  86. 0x405c,0x33f2,0x8a58,0x1b86,
  87. 0x40a1,0x76fa,0x0e55,0x35fa,
  88. 0x40b2,0xec10,0x2442,0x040c
  89. };
  90. #endif
  91. #ifdef ANSIPROT
  92. extern double fabs ( double );
  93. extern double exp ( double );
  94. extern double polevl ( double, void *, int );
  95. extern double p1evl ( double, void *, int );
  96. #else
  97. double fabs(), exp(), polevl(), p1evl();
  98. #endif
  99. extern double MAXLOG;
  100. double tanh(x)
  101. double x;
  102. {
  103. double s, z;
  104. #ifdef MINUSZERO
  105. if( x == 0.0 )
  106. return(x);
  107. #endif
  108. z = fabs(x);
  109. if( z > 0.5 * MAXLOG )
  110. {
  111. if( x > 0 )
  112. return( 1.0 );
  113. else
  114. return( -1.0 );
  115. }
  116. if( z >= 0.625 )
  117. {
  118. s = exp(2.0*z);
  119. z = 1.0 - 2.0/(s + 1.0);
  120. if( x < 0 )
  121. z = -z;
  122. }
  123. else
  124. {
  125. if( x == 0.0 )
  126. return(x);
  127. s = x * x;
  128. z = polevl( s, P, 2 )/p1evl(s, Q, 3);
  129. z = x * s * z;
  130. z = x + z;
  131. }
  132. return( z );
  133. }