acosh.c 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167
  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 <math.h>
  50. #ifdef UNK
  51. static double P[] = {
  52. 1.18801130533544501356E2,
  53. 3.94726656571334401102E3,
  54. 3.43989375926195455866E4,
  55. 1.08102874834699867335E5,
  56. 1.10855947270161294369E5
  57. };
  58. static double Q[] = {
  59. /* 1.00000000000000000000E0,*/
  60. 1.86145380837903397292E2,
  61. 4.15352677227719831579E3,
  62. 2.97683430363289370382E4,
  63. 8.29725251988426222434E4,
  64. 7.83869920495893927727E4
  65. };
  66. #endif
  67. #ifdef DEC
  68. static unsigned short P[] = {
  69. 0041755,0115055,0144002,0146444,
  70. 0043166,0132103,0155150,0150302,
  71. 0044006,0057360,0003021,0162753,
  72. 0044323,0021557,0175225,0056253,
  73. 0044330,0101771,0040046,0006636
  74. };
  75. static unsigned short Q[] = {
  76. /*0040200,0000000,0000000,0000000,*/
  77. 0042072,0022467,0126670,0041232,
  78. 0043201,0146066,0152142,0034015,
  79. 0043750,0110257,0121165,0026100,
  80. 0044242,0007103,0034667,0033173,
  81. 0044231,0014576,0175573,0017472
  82. };
  83. #endif
  84. #ifdef IBMPC
  85. static unsigned short P[] = {
  86. 0x59a4,0xb900,0xb345,0x405d,
  87. 0x1a18,0x7b4d,0xd688,0x40ae,
  88. 0x3cbd,0x00c2,0xcbde,0x40e0,
  89. 0xab95,0xff52,0x646d,0x40fa,
  90. 0xc1b4,0x2804,0x107f,0x40fb
  91. };
  92. static unsigned short Q[] = {
  93. /*0x0000,0x0000,0x0000,0x3ff0,*/
  94. 0x0853,0xf5b7,0x44a6,0x4067,
  95. 0x4702,0xda8c,0x3986,0x40b0,
  96. 0xa588,0xf44e,0x1215,0x40dd,
  97. 0xe6cf,0x6736,0x41c8,0x40f4,
  98. 0x63e7,0xdf6f,0x232f,0x40f3
  99. };
  100. #endif
  101. #ifdef MIEEE
  102. static unsigned short P[] = {
  103. 0x405d,0xb345,0xb900,0x59a4,
  104. 0x40ae,0xd688,0x7b4d,0x1a18,
  105. 0x40e0,0xcbde,0x00c2,0x3cbd,
  106. 0x40fa,0x646d,0xff52,0xab95,
  107. 0x40fb,0x107f,0x2804,0xc1b4
  108. };
  109. static unsigned short Q[] = {
  110. 0x4067,0x44a6,0xf5b7,0x0853,
  111. 0x40b0,0x3986,0xda8c,0x4702,
  112. 0x40dd,0x1215,0xf44e,0xa588,
  113. 0x40f4,0x41c8,0x6736,0xe6cf,
  114. 0x40f3,0x232f,0xdf6f,0x63e7,
  115. };
  116. #endif
  117. #ifdef ANSIPROT
  118. extern double polevl ( double, void *, int );
  119. extern double p1evl ( double, void *, int );
  120. extern double log ( double );
  121. extern double sqrt ( double );
  122. #else
  123. double log(), sqrt(), polevl(), p1evl();
  124. #endif
  125. extern double LOGE2, INFINITY, NAN;
  126. double acosh(x)
  127. double x;
  128. {
  129. double a, z;
  130. if( x < 1.0 )
  131. {
  132. mtherr( "acosh", DOMAIN );
  133. return(NAN);
  134. }
  135. if( x > 1.0e8 )
  136. {
  137. #ifdef INFINITIES
  138. if( x == INFINITY )
  139. return( INFINITY );
  140. #endif
  141. return( log(x) + LOGE2 );
  142. }
  143. z = x - 1.0;
  144. if( z < 0.5 )
  145. {
  146. a = sqrt(z) * (polevl(z, P, 4) / p1evl(z, Q, 5) );
  147. return( a );
  148. }
  149. a = sqrt( z*(x+1.0) );
  150. return( log(x + a) );
  151. }