sinh.c 1.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081
  1. /*
  2. * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  3. * See the copyright notice in the ACK home directory, in the file "Copyright".
  4. *
  5. * Author: Ceriel J.H. Jacobs
  6. */
  7. /* $Id$ */
  8. #include <math.h>
  9. #include <float.h>
  10. #include <errno.h>
  11. #include "localmath.h"
  12. static double
  13. sinh_cosh(double x, int cosh_flag)
  14. {
  15. /* Algorithm and coefficients from:
  16. "Software manual for the elementary functions"
  17. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  18. */
  19. static double p[] = {
  20. -0.35181283430177117881e+6,
  21. -0.11563521196851768270e+5,
  22. -0.16375798202630751372e+3,
  23. -0.78966127417357099479e+0
  24. };
  25. static double q[] = {
  26. -0.21108770058106271242e+7,
  27. 0.36162723109421836460e+5,
  28. -0.27773523119650701167e+3,
  29. 1.0
  30. };
  31. int negative = x < 0;
  32. double y = negative ? -x : x;
  33. if (__IsNan(x)) {
  34. errno = EDOM;
  35. return x;
  36. }
  37. if (! cosh_flag && y <= 1.0) {
  38. /* ??? check for underflow ??? */
  39. y = y * y;
  40. return x + x * y * POLYNOM3(y, p)/POLYNOM3(y,q);
  41. }
  42. if (y >= M_LN_MAX_D) {
  43. /* exp(y) would cause overflow */
  44. #define LNV 0.69316101074218750000e+0
  45. #define VD2M1 0.52820835025874852469e-4
  46. double w = y - LNV;
  47. if (w < M_LN_MAX_D+M_LN2-LNV) {
  48. x = exp(w);
  49. x += VD2M1 * x;
  50. }
  51. else {
  52. errno = ERANGE;
  53. x = HUGE_VAL;
  54. }
  55. }
  56. else {
  57. double z = exp(y);
  58. x = 0.5 * (z + (cosh_flag ? 1.0 : -1.0)/z);
  59. }
  60. return negative ? -x : x;
  61. }
  62. double
  63. sinh(double x)
  64. {
  65. return sinh_cosh(x, 0);
  66. }
  67. double
  68. cosh(double x)
  69. {
  70. if (x < 0) x = -x;
  71. return sinh_cosh(x, 1);
  72. }