acosh.c 3.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138
  1. /* acosh.c
  2. *
  3. * Inverse hyperbolic cosine
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, acosh();
  10. *
  11. * y = acosh( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns inverse hyperbolic cosine of argument.
  18. *
  19. * If 1 <= x < 1.5, a rational approximation
  20. *
  21. * sqrt(z) * P(z)/Q(z)
  22. *
  23. * where z = x-1, is used. Otherwise,
  24. *
  25. * acosh(x) = log( x + sqrt( (x-1)(x+1) ).
  26. *
  27. *
  28. *
  29. * ACCURACY:
  30. *
  31. * Relative error:
  32. * arithmetic domain # trials peak rms
  33. * DEC 1,3 30000 4.2e-17 1.1e-17
  34. * IEEE 1,3 30000 4.6e-16 8.7e-17
  35. *
  36. *
  37. * ERROR MESSAGES:
  38. *
  39. * message condition value returned
  40. * acosh domain |x| < 1 NAN
  41. *
  42. */
  43. /* acosh.c */
  44. /*
  45. Cephes Math Library Release 2.8: June, 2000
  46. Copyright 1984, 1995, 2000 by Stephen L. Moshier
  47. */
  48. /* acosh(z) = sqrt(x) * R(x), z = x + 1, interval 0 < x < 0.5 */
  49. #include "mconf.h"
  50. #ifdef UNK
  51. static double P[] = {1.18801130533544501356E2, 3.94726656571334401102E3,
  52. 3.43989375926195455866E4, 1.08102874834699867335E5,
  53. 1.10855947270161294369E5};
  54. static double Q[] = {
  55. /* 1.00000000000000000000E0,*/
  56. 1.86145380837903397292E2, 4.15352677227719831579E3,
  57. 2.97683430363289370382E4, 8.29725251988426222434E4,
  58. 7.83869920495893927727E4};
  59. #endif
  60. #ifdef DEC
  61. static unsigned short P[] = {0041755, 0115055, 0144002, 0146444, 0043166,
  62. 0132103, 0155150, 0150302, 0044006, 0057360,
  63. 0003021, 0162753, 0044323, 0021557, 0175225,
  64. 0056253, 0044330, 0101771, 0040046, 0006636};
  65. static unsigned short Q[] = {
  66. /*0040200,0000000,0000000,0000000,*/
  67. 0042072, 0022467, 0126670, 0041232, 0043201, 0146066, 0152142,
  68. 0034015, 0043750, 0110257, 0121165, 0026100, 0044242, 0007103,
  69. 0034667, 0033173, 0044231, 0014576, 0175573, 0017472};
  70. #endif
  71. #ifdef IBMPC
  72. static unsigned short P[] = {0x59a4, 0xb900, 0xb345, 0x405d, 0x1a18,
  73. 0x7b4d, 0xd688, 0x40ae, 0x3cbd, 0x00c2,
  74. 0xcbde, 0x40e0, 0xab95, 0xff52, 0x646d,
  75. 0x40fa, 0xc1b4, 0x2804, 0x107f, 0x40fb};
  76. static unsigned short Q[] = {
  77. /*0x0000,0x0000,0x0000,0x3ff0,*/
  78. 0x0853, 0xf5b7, 0x44a6, 0x4067, 0x4702, 0xda8c, 0x3986,
  79. 0x40b0, 0xa588, 0xf44e, 0x1215, 0x40dd, 0xe6cf, 0x6736,
  80. 0x41c8, 0x40f4, 0x63e7, 0xdf6f, 0x232f, 0x40f3};
  81. #endif
  82. #ifdef MIEEE
  83. static unsigned short P[] = {0x405d, 0xb345, 0xb900, 0x59a4, 0x40ae,
  84. 0xd688, 0x7b4d, 0x1a18, 0x40e0, 0xcbde,
  85. 0x00c2, 0x3cbd, 0x40fa, 0x646d, 0xff52,
  86. 0xab95, 0x40fb, 0x107f, 0x2804, 0xc1b4};
  87. static unsigned short Q[] = {
  88. 0x4067, 0x44a6, 0xf5b7, 0x0853, 0x40b0, 0x3986, 0xda8c,
  89. 0x4702, 0x40dd, 0x1215, 0xf44e, 0xa588, 0x40f4, 0x41c8,
  90. 0x6736, 0xe6cf, 0x40f3, 0x232f, 0xdf6f, 0x63e7,
  91. };
  92. #endif
  93. #ifdef ANSIPROT
  94. extern double polevl(double, void *, int);
  95. extern double p1evl(double, void *, int);
  96. extern double log(double);
  97. extern double sqrt(double);
  98. #else
  99. double log(), sqrt(), polevl(), p1evl();
  100. #endif
  101. extern double LOGE2, INFINITY, NAN;
  102. double acosh(x) double x;
  103. {
  104. double a, z;
  105. if (x < 1.0) {
  106. mtherr("acosh", DOMAIN);
  107. return (NAN);
  108. }
  109. if (x > 1.0e8) {
  110. #ifdef INFINITIES
  111. if (x == INFINITY)
  112. return (INFINITY);
  113. #endif
  114. return (log(x) + LOGE2);
  115. }
  116. z = x - 1.0;
  117. if (z < 0.5) {
  118. a = sqrt(z) * (polevl(z, P, 4) / p1evl(z, Q, 5));
  119. return (a);
  120. }
  121. a = sqrt(z * (x + 1.0));
  122. return (log(x + a));
  123. }