asinh.c 3.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133
  1. /* asinh.c
  2. *
  3. * Inverse hyperbolic sine
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, asinh();
  10. *
  11. * y = asinh( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns inverse hyperbolic sine of argument.
  18. *
  19. * If |x| < 0.5, the function is approximated by a rational
  20. * form x + x**3 P(x)/Q(x). Otherwise,
  21. *
  22. * asinh(x) = log( x + sqrt(1 + x*x) ).
  23. *
  24. *
  25. *
  26. * ACCURACY:
  27. *
  28. * Relative error:
  29. * arithmetic domain # trials peak rms
  30. * DEC -3,3 75000 4.6e-17 1.1e-17
  31. * IEEE -1,1 30000 3.7e-16 7.8e-17
  32. * IEEE 1,3 30000 2.5e-16 6.7e-17
  33. *
  34. */
  35. /* asinh.c */
  36. /*
  37. Cephes Math Library Release 2.8: June, 2000
  38. Copyright 1984, 1995, 2000 by Stephen L. Moshier
  39. */
  40. #include "mconf.h"
  41. #ifdef UNK
  42. static double P[] = {-4.33231683752342103572E-3, -5.91750212056387121207E-1,
  43. -4.37390226194356683570E0, -9.09030533308377316566E0,
  44. -5.56682227230859640450E0};
  45. static double Q[] = {
  46. /* 1.00000000000000000000E0,*/
  47. 1.28757002067426453537E1, 4.86042483805291788324E1,
  48. 6.95722521337257608734E1, 3.34009336338516356383E1};
  49. #endif
  50. #ifdef DEC
  51. static unsigned short P[] = {0136215, 0173033, 0110410, 0105475, 0140027,
  52. 0076361, 0020056, 0164520, 0140613, 0173401,
  53. 0160136, 0053142, 0141021, 0070744, 0000503,
  54. 0176261, 0140662, 0021550, 0073106, 0133351};
  55. static unsigned short Q[] = {
  56. /* 0040200,0000000,0000000,0000000,*/
  57. 0041116, 0001336, 0034120, 0173054, 0041502, 0065300, 0013144, 0021231,
  58. 0041613, 0022376, 0035516, 0153063, 0041405, 0115216, 0054265, 0004557};
  59. #endif
  60. #ifdef IBMPC
  61. static unsigned short P[] = {0x1168, 0x7221, 0xbec3, 0xbf71, 0xdd2a,
  62. 0x2405, 0xef9e, 0xbfe2, 0xcacc, 0x3c0b,
  63. 0x7ee0, 0xc011, 0x7f96, 0x8028, 0x2e3c,
  64. 0xc022, 0xd6dd, 0x0ec8, 0x446d, 0xc016};
  65. static unsigned short Q[] = {
  66. /* 0x0000,0x0000,0x0000,0x3ff0,*/
  67. 0x1ec5, 0xc70a, 0xc05b, 0x4029, 0x8453, 0x02cc, 0x4d58, 0x4048,
  68. 0xdac6, 0xc769, 0x649f, 0x4051, 0xa12e, 0xcb16, 0xb351, 0x4040};
  69. #endif
  70. #ifdef MIEEE
  71. static unsigned short P[] = {0xbf71, 0xbec3, 0x7221, 0x1168, 0xbfe2,
  72. 0xef9e, 0x2405, 0xdd2a, 0xc011, 0x7ee0,
  73. 0x3c0b, 0xcacc, 0xc022, 0x2e3c, 0x8028,
  74. 0x7f96, 0xc016, 0x446d, 0x0ec8, 0xd6dd};
  75. static unsigned short Q[] = {0x4029, 0xc05b, 0xc70a, 0x1ec5, 0x4048, 0x4d58,
  76. 0x02cc, 0x8453, 0x4051, 0x649f, 0xc769, 0xdac6,
  77. 0x4040, 0xb351, 0xcb16, 0xa12e};
  78. #endif
  79. #ifdef ANSIPROT
  80. extern double polevl(double, void *, int);
  81. extern double p1evl(double, void *, int);
  82. extern double sqrt(double);
  83. extern double log(double);
  84. #else
  85. double log(), sqrt(), polevl(), p1evl();
  86. #endif
  87. extern double LOGE2, INFINITY;
  88. double asinh(xx) double xx;
  89. {
  90. double a, z, x;
  91. int sign;
  92. #ifdef MINUSZERO
  93. if (xx == 0.0)
  94. return (xx);
  95. #endif
  96. if (xx < 0.0) {
  97. sign = -1;
  98. x = -xx;
  99. } else {
  100. sign = 1;
  101. x = xx;
  102. }
  103. if (x > 1.0e8) {
  104. #ifdef INFINITIES
  105. if (x == INFINITY)
  106. return (xx);
  107. #endif
  108. return (sign * (log(x) + LOGE2));
  109. }
  110. z = x * x;
  111. if (x < 0.5) {
  112. a = (polevl(z, P, 4) / p1evl(z, Q, 4)) * z;
  113. a = a * x + x;
  114. if (sign < 0)
  115. a = -a;
  116. return (a);
  117. }
  118. a = sqrt(z + 1.0);
  119. return (sign * log(x + a));
  120. }