Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /*                                                      tanhl.c
  2.  *
  3.  *      Hyperbolic tangent, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, tanhl();
  10.  *
  11.  * y = tanhl( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns hyperbolic tangent of argument in the range MINLOGL to
  18.  * MAXLOGL.
  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.  *    IEEE      -2,2        30000       1.3e-19     2.4e-20
  32.  *
  33.  */
  35. /*
  36. Cephes Math Library Release 2.7:  May, 1998
  37. Copyright 1984, 1987, 1989, 1998 by Stephen L. Moshier
  38. */
  39.  
  40. /*
  41. Modified for mingw
  42. 2002-07-22 Danny Smith <dannysmith@users.sourceforge.net>
  43. */
  44.  
  45. #ifdef __MINGW32__
  46. #include "cephes_mconf.h"
  47. #else
  48. #include "mconf.h"
  49. #endif
  50.  
  51. #ifndef _SET_ERRNO
  52. #define _SET_ERRNO(x)
  53. #endif
  54.  
  55. #ifdef UNK
  56. static long double P[] = {
  57. -6.8473739392677100872869E-5L,
  58. -9.5658283111794641589011E-1L,
  59. -8.4053568599672284488465E1L,
  60. -1.3080425704712825945553E3L,
  61. };
  62. static long double Q[] = {
  63. /* 1.0000000000000000000000E0L,*/
  64.  9.6259501838840336946872E1L,
  65.  1.8218117903645559060232E3L,
  66.  3.9241277114138477845780E3L,
  67. };
  68. #endif
  69.  
  70. #ifdef IBMPC
  71. static unsigned short P[] = {
  72. 0xd2a4,0x1b0c,0x8f15,0x8f99,0xbff1, XPD
  73. 0x5959,0x9111,0x9cc7,0xf4e2,0xbffe, XPD
  74. 0xb576,0xef5e,0x6d57,0xa81b,0xc005, XPD
  75. 0xe3be,0xbfbd,0x5cbc,0xa381,0xc009, XPD
  76. };
  77. static unsigned short Q[] = {
  78. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  79. 0x687f,0xce24,0xdd6c,0xc084,0x4005, XPD
  80. 0x3793,0xc95f,0xfa2f,0xe3b9,0x4009, XPD
  81. 0xd5a2,0x1f9c,0x0b1b,0xf542,0x400a, XPD
  82. };
  83. #endif
  84.  
  85. #ifdef MIEEE
  86. static long P[] = {
  87. 0xbff10000,0x8f998f15,0x1b0cd2a4,
  88. 0xbffe0000,0xf4e29cc7,0x91115959,
  89. 0xc0050000,0xa81b6d57,0xef5eb576,
  90. 0xc0090000,0xa3815cbc,0xbfbde3be,
  91. };
  92. static long Q[] = {
  93. /*0x3fff0000,0x80000000,0x00000000,*/
  94. 0x40050000,0xc084dd6c,0xce24687f,
  95. 0x40090000,0xe3b9fa2f,0xc95f3793,
  96. 0x400a0000,0xf5420b1b,0x1f9cd5a2,
  97. };
  98. #endif
  99.  
  100. #ifndef __MINGW32__
  101. extern long double MAXLOGL;
  102. #ifdef ANSIPROT
  103. extern long double fabsl ( long double );
  104. extern long double expl ( long double );
  105. extern long double polevll ( long double, void *, int );
  106. extern long double p1evll ( long double, void *, int );
  107. #else
  108. long double fabsl(), expl(), polevll(), p1evll();
  109. #endif
  110. #endif /* __MINGW32__ */
  111.  
  112. long double tanhl(x)
  113. long double x;
  114. {
  115. long double s, z;
  116.  
  117. #ifdef MINUSZERO
  118. if( x == 0.0L )
  119.         return(x);
  120. #endif
  121. if (isnanl(x))
  122.         {
  123.         _SET_ERRNO (EDOM);
  124.         return x;
  125.         }
  126.  
  127. z = fabsl(x);
  128. if( z > 0.5L * MAXLOGL )
  129.         {
  130.         _SET_ERRNO (ERANGE);
  131.         if( x > 0 )
  132.                 return( 1.0L );
  133.         else
  134.                 return( -1.0L );
  135.         }
  136. if( z >= 0.625L )
  137.         {
  138.         s = expl(2.0*z);
  139.         z =  1.0L  - 2.0/(s + 1.0L);
  140.         if( x < 0 )
  141.                 z = -z;
  142.         }
  143. else
  144.         {
  145.         s = x * x;
  146.         z = polevll( s, P, 3 )/p1evll(s, Q, 3);
  147.         z = x * s * z;
  148.         z = x + z;
  149.         }
  150. return( z );
  151. }
  152.