sinh.c 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122
  1. /* sinh.c
  2. *
  3. * Hyperbolic sine
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, sinh();
  10. *
  11. * y = sinh( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns hyperbolic sine of argument in the range MINLOG to
  18. * MAXLOG.
  19. *
  20. * The range is partitioned into two segments. If |x| <= 1, a
  21. * rational function of the form x + x**3 P(x)/Q(x) is employed.
  22. * Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2.
  23. *
  24. *
  25. *
  26. * ACCURACY:
  27. *
  28. * Relative error:
  29. * arithmetic domain # trials peak rms
  30. * DEC +- 88 50000 4.0e-17 7.7e-18
  31. * IEEE +-MAXLOG 30000 2.6e-16 5.7e-17
  32. *
  33. */
  34. /*
  35. Cephes Math Library Release 2.8: June, 2000
  36. Copyright 1984, 1995, 2000 by Stephen L. Moshier
  37. */
  38. #include "mconf.h"
  39. #ifdef UNK
  40. static double P[] = {-7.89474443963537015605E-1, -1.63725857525983828727E2,
  41. -1.15614435765005216044E4, -3.51754964808151394800E5};
  42. static double Q[] = {
  43. /* 1.00000000000000000000E0,*/
  44. -2.77711081420602794433E2, 3.61578279834431989373E4,
  45. -2.11052978884890840399E6};
  46. #endif
  47. #ifdef DEC
  48. static unsigned short P[] = {
  49. 0140112, 0015377, 0042731, 0163255, 0142043, 0134721, 0146177, 0123761,
  50. 0143464, 0122706, 0034353, 0006017, 0144653, 0140536, 0157665, 0054045};
  51. static unsigned short Q[] = {
  52. /*0040200,0000000,0000000,0000000,*/
  53. 0142212, 0155404, 0133513, 0022040, 0044015, 0036723,
  54. 0173271, 0011053, 0145400, 0150407, 0023710, 0001034};
  55. #endif
  56. #ifdef IBMPC
  57. static unsigned short P[] = {0x3cd6, 0xe8bb, 0x435f, 0xbfe9, 0xf4fe, 0x398f,
  58. 0x773a, 0xc064, 0x6182, 0xc71d, 0x94b8, 0xc0c6,
  59. 0xab05, 0xdbf6, 0x782b, 0xc115};
  60. static unsigned short Q[] = {
  61. /*0x0000,0x0000,0x0000,0x3ff0,*/
  62. 0x6484, 0x96e9, 0x5b60, 0xc071, 0x2245, 0x7ed7,
  63. 0xa7ba, 0x40e1, 0x0044, 0xe4f9, 0x1a20, 0xc140};
  64. #endif
  65. #ifdef MIEEE
  66. static unsigned short P[] = {0xbfe9, 0x435f, 0xe8bb, 0x3cd6, 0xc064, 0x773a,
  67. 0x398f, 0xf4fe, 0xc0c6, 0x94b8, 0xc71d, 0x6182,
  68. 0xc115, 0x782b, 0xdbf6, 0xab05};
  69. static unsigned short Q[] = {0xc071, 0x5b60, 0x96e9, 0x6484, 0x40e1, 0xa7ba,
  70. 0x7ed7, 0x2245, 0xc140, 0x1a20, 0xe4f9, 0x0044};
  71. #endif
  72. #ifdef ANSIPROT
  73. extern double fabs(double);
  74. extern double exp(double);
  75. extern double polevl(double, void *, int);
  76. extern double p1evl(double, void *, int);
  77. #else
  78. double fabs(), exp(), polevl(), p1evl();
  79. #endif
  80. extern double INFINITY, MINLOG, MAXLOG, LOGE2;
  81. double sinh(x) double x;
  82. {
  83. double a;
  84. #ifdef MINUSZERO
  85. if (x == 0.0)
  86. return (x);
  87. #endif
  88. a = fabs(x);
  89. if ((x > (MAXLOG + LOGE2)) || (x > -(MINLOG - LOGE2))) {
  90. mtherr("sinh", DOMAIN);
  91. if (x > 0)
  92. return (INFINITY);
  93. else
  94. return (-INFINITY);
  95. }
  96. if (a > 1.0) {
  97. if (a >= (MAXLOG - LOGE2)) {
  98. a = exp(0.5 * a);
  99. a = (0.5 * a) * a;
  100. if (x < 0)
  101. a = -a;
  102. return (a);
  103. }
  104. a = exp(a);
  105. a = 0.5 * a - (0.5 / a);
  106. if (x < 0)
  107. a = -a;
  108. return (a);
  109. }
  110. a *= a;
  111. return (x + x * a * (polevl(a, P, 3) / p1evl(a, Q, 3)));
  112. }