exp2.c 3.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171
  1. /* exp2.c
  2. *
  3. * Base 2 exponential function
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, exp2();
  10. *
  11. * y = exp2( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns 2 raised to the x power.
  18. *
  19. * Range reduction is accomplished by separating the argument
  20. * into an integer k and fraction f such that
  21. * x k f
  22. * 2 = 2 2.
  23. *
  24. * A Pade' form
  25. *
  26. * 1 + 2x P(x**2) / (Q(x**2) - x P(x**2) )
  27. *
  28. * approximates 2**x in the basic range [-0.5, 0.5].
  29. *
  30. *
  31. * ACCURACY:
  32. *
  33. * Relative error:
  34. * arithmetic domain # trials peak rms
  35. * IEEE -1022,+1024 30000 1.8e-16 5.4e-17
  36. *
  37. *
  38. * See exp.c for comments on error amplification.
  39. *
  40. *
  41. * ERROR MESSAGES:
  42. *
  43. * message condition value returned
  44. * exp underflow x < -MAXL2 0.0
  45. * exp overflow x > MAXL2 MAXNUM
  46. *
  47. * For DEC arithmetic, MAXL2 = 127.
  48. * For IEEE arithmetic, MAXL2 = 1024.
  49. */
  50. /*
  51. Cephes Math Library Release 2.8: June, 2000
  52. Copyright 1984, 1995, 2000 by Stephen L. Moshier
  53. */
  54. #include "mconf.h"
  55. #ifdef UNK
  56. static double P[] = {
  57. 2.30933477057345225087E-2,
  58. 2.02020656693165307700E1,
  59. 1.51390680115615096133E3,
  60. };
  61. static double Q[] = {
  62. /* 1.00000000000000000000E0,*/
  63. 2.33184211722314911771E2,
  64. 4.36821166879210612817E3,
  65. };
  66. #define MAXL2 1024.0
  67. #define MINL2 -1024.0
  68. #endif
  69. #ifdef DEC
  70. static unsigned short P[] = {
  71. 0036675, 0027102, 0122327, 0053227, 0041241, 0116724,
  72. 0115412, 0157355, 0042675, 0036404, 0101733, 0132226,
  73. };
  74. static unsigned short Q[] = {
  75. /*0040200,0000000,0000000,0000000,*/
  76. 0042151, 0027450, 0077732, 0160744, 0043210, 0100661, 0077550, 0056560,
  77. };
  78. #define MAXL2 127.0
  79. #define MINL2 -127.0
  80. #endif
  81. #ifdef IBMPC
  82. static unsigned short P[] = {
  83. 0xead3, 0x549a, 0xa5c8, 0x3f97, 0x5bde, 0x9361,
  84. 0x33ba, 0x4034, 0x7693, 0x907b, 0xa7a0, 0x4097,
  85. };
  86. static unsigned short Q[] = {
  87. /*0x0000,0x0000,0x0000,0x3ff0,*/
  88. 0x5c3c, 0x0ffb, 0x25e5, 0x406d, 0x0bae, 0x2fed, 0x1036, 0x40b1,
  89. };
  90. #define MAXL2 1024.0
  91. #define MINL2 -1022.0
  92. #endif
  93. #ifdef MIEEE
  94. static unsigned short P[] = {
  95. 0x3f97, 0xa5c8, 0x549a, 0xead3, 0x4034, 0x33ba,
  96. 0x9361, 0x5bde, 0x4097, 0xa7a0, 0x907b, 0x7693,
  97. };
  98. static unsigned short Q[] = {
  99. /*0x3ff0,0x0000,0x0000,0x0000,*/
  100. 0x406d, 0x25e5, 0x0ffb, 0x5c3c, 0x40b1, 0x1036, 0x2fed, 0x0bae,
  101. };
  102. #define MAXL2 1024.0
  103. #define MINL2 -1022.0
  104. #endif
  105. #ifdef ANSIPROT
  106. extern double polevl(double, void *, int);
  107. extern double p1evl(double, void *, int);
  108. extern double floor(double);
  109. extern double ldexp(double, int);
  110. extern int isnan(double);
  111. extern int isfinite(double);
  112. #else
  113. double polevl(), p1evl(), floor(), ldexp();
  114. int isnan(), isfinite();
  115. #endif
  116. #ifdef INFINITIES
  117. extern double INFINITY;
  118. #endif
  119. extern double MAXNUM;
  120. double exp2(x) double x;
  121. {
  122. double px, xx;
  123. short n;
  124. #ifdef NANS
  125. if (isnan(x))
  126. return (x);
  127. #endif
  128. if (x > MAXL2) {
  129. #ifdef INFINITIES
  130. return (INFINITY);
  131. #else
  132. mtherr("exp2", OVERFLOW);
  133. return (MAXNUM);
  134. #endif
  135. }
  136. if (x < MINL2) {
  137. #ifndef INFINITIES
  138. mtherr("exp2", UNDERFLOW);
  139. #endif
  140. return (0.0);
  141. }
  142. xx = x; /* save x */
  143. /* separate into integer and fractional parts */
  144. px = floor(x + 0.5);
  145. n = px;
  146. x = x - px;
  147. /* rational approximation
  148. * exp2(x) = 1 + 2xP(xx)/(Q(xx) - P(xx))
  149. * where xx = x**2
  150. */
  151. xx = x * x;
  152. px = x * polevl(xx, P, 2);
  153. x = px / (p1evl(xx, Q, 2) - px);
  154. x = 1.0 + ldexp(x, 1);
  155. /* scale by power of 2 */
  156. x = ldexp(x, n);
  157. return (x);
  158. }