exp10.c 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246
  1. /* exp10.c
  2. *
  3. * Base 10 exponential function
  4. * (Common antilogarithm)
  5. *
  6. *
  7. *
  8. * SYNOPSIS:
  9. *
  10. * double x, y, exp10();
  11. *
  12. * y = exp10( x );
  13. *
  14. *
  15. *
  16. * DESCRIPTION:
  17. *
  18. * Returns 10 raised to the x power.
  19. *
  20. * Range reduction is accomplished by expressing the argument
  21. * as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
  22. * The Pade' form
  23. *
  24. * 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
  25. *
  26. * is used to approximate 10**f.
  27. *
  28. *
  29. *
  30. * ACCURACY:
  31. *
  32. * Relative error:
  33. * arithmetic domain # trials peak rms
  34. * IEEE -307,+307 30000 2.2e-16 5.5e-17
  35. * Test result from an earlier version (2.1):
  36. * DEC -38,+38 70000 3.1e-17 7.0e-18
  37. *
  38. * ERROR MESSAGES:
  39. *
  40. * message condition value returned
  41. * exp10 underflow x < -MAXL10 0.0
  42. * exp10 overflow x > MAXL10 MAXNUM
  43. *
  44. * DEC arithmetic: MAXL10 = 38.230809449325611792.
  45. * IEEE arithmetic: MAXL10 = 308.2547155599167.
  46. *
  47. */
  48. /*
  49. Cephes Math Library Release 2.8: June, 2000
  50. Copyright 1984, 1991, 2000 by Stephen L. Moshier
  51. */
  52. #include "mconf.h"
  53. #ifdef UNK
  54. static double P[] = {
  55. 4.09962519798587023075E-2,
  56. 1.17452732554344059015E1,
  57. 4.06717289936872725516E2,
  58. 2.39423741207388267439E3,
  59. };
  60. static double Q[] = {
  61. /* 1.00000000000000000000E0,*/
  62. 8.50936160849306532625E1,
  63. 1.27209271178345121210E3,
  64. 2.07960819286001865907E3,
  65. };
  66. /* static double LOG102 = 3.01029995663981195214e-1; */
  67. static double LOG210 = 3.32192809488736234787e0;
  68. static double LG102A = 3.01025390625000000000E-1;
  69. static double LG102B = 4.60503898119521373889E-6;
  70. /* static double MAXL10 = 38.230809449325611792; */
  71. static double MAXL10 = 308.2547155599167;
  72. #endif
  73. #ifdef DEC
  74. static unsigned short P[] = {
  75. 0037047, 0165657, 0114061, 0067234, 0041073, 0166243, 0123052, 0144643,
  76. 0042313, 0055720, 0024032, 0047443, 0043025, 0121714, 0070232, 0050007,
  77. };
  78. static unsigned short Q[] = {
  79. /*0040200,0000000,0000000,0000000,*/
  80. 0041652, 0027756, 0071216, 0050075, 0042637, 0001367,
  81. 0077263, 0136017, 0043001, 0174673, 0024157, 0133416,
  82. };
  83. /*
  84. static unsigned short L102[] = {0037632,0020232,0102373,0147770};
  85. #define LOG102 *(double *)L102
  86. */
  87. static unsigned short L210[] = {0040524, 0115170, 0045715, 0015613};
  88. #define LOG210 *(double *)L210
  89. static unsigned short L102A[] = {
  90. 0037632,
  91. 0020000,
  92. 0000000,
  93. 0000000,
  94. };
  95. #define LG102A *(double *)L102A
  96. static unsigned short L102B[] = {
  97. 0033632,
  98. 0102373,
  99. 0147767,
  100. 0114220,
  101. };
  102. #define LG102B *(double *)L102B
  103. static unsigned short MXL[] = {
  104. 0041430,
  105. 0166131,
  106. 0047761,
  107. 0154130,
  108. };
  109. #define MAXL10 (*(double *)MXL)
  110. #endif
  111. #ifdef IBMPC
  112. static unsigned short P[] = {
  113. 0x2dd4, 0xf306, 0xfd75, 0x3fa4, 0x5934, 0x74c5, 0x7d94, 0x4027,
  114. 0x49e4, 0x0503, 0x6b7a, 0x4079, 0x4a01, 0x8e13, 0xb479, 0x40a2,
  115. };
  116. static unsigned short Q[] = {
  117. /*0x0000,0x0000,0x0000,0x3ff0,*/
  118. 0xca08, 0xce51, 0x45fd, 0x4055, 0x7782, 0xefd6,
  119. 0xe05e, 0x4093, 0xf6e2, 0x650d, 0x3f37, 0x40a0,
  120. };
  121. /*
  122. static unsigned short L102[] = {0x79ff,0x509f,0x4413,0x3fd3};
  123. #define LOG102 *(double *)L102
  124. */
  125. static unsigned short L210[] = {0xa371, 0x0979, 0x934f, 0x400a};
  126. #define LOG210 *(double *)L210
  127. static unsigned short L102A[] = {
  128. 0x0000,
  129. 0x0000,
  130. 0x4400,
  131. 0x3fd3,
  132. };
  133. #define LG102A *(double *)L102A
  134. static unsigned short L102B[] = {
  135. 0xf312,
  136. 0x79fe,
  137. 0x509f,
  138. 0x3ed3,
  139. };
  140. #define LG102B *(double *)L102B
  141. static double MAXL10 = 308.2547155599167;
  142. #endif
  143. #ifdef MIEEE
  144. static unsigned short P[] = {
  145. 0x3fa4, 0xfd75, 0xf306, 0x2dd4, 0x4027, 0x7d94, 0x74c5, 0x5934,
  146. 0x4079, 0x6b7a, 0x0503, 0x49e4, 0x40a2, 0xb479, 0x8e13, 0x4a01,
  147. };
  148. static unsigned short Q[] = {
  149. /*0x3ff0,0x0000,0x0000,0x0000,*/
  150. 0x4055, 0x45fd, 0xce51, 0xca08, 0x4093, 0xe05e,
  151. 0xefd6, 0x7782, 0x40a0, 0x3f37, 0x650d, 0xf6e2,
  152. };
  153. /*
  154. static unsigned short L102[] = {0x3fd3,0x4413,0x509f,0x79ff};
  155. #define LOG102 *(double *)L102
  156. */
  157. static unsigned short L210[] = {0x400a, 0x934f, 0x0979, 0xa371};
  158. #define LOG210 *(double *)L210
  159. static unsigned short L102A[] = {
  160. 0x3fd3,
  161. 0x4400,
  162. 0x0000,
  163. 0x0000,
  164. };
  165. #define LG102A *(double *)L102A
  166. static unsigned short L102B[] = {
  167. 0x3ed3,
  168. 0x509f,
  169. 0x79fe,
  170. 0xf312,
  171. };
  172. #define LG102B *(double *)L102B
  173. static double MAXL10 = 308.2547155599167;
  174. #endif
  175. #ifdef ANSIPROT
  176. extern double floor(double);
  177. extern double ldexp(double, int);
  178. extern double polevl(double, void *, int);
  179. extern double p1evl(double, void *, int);
  180. extern int isnan(double);
  181. extern int isfinite(double);
  182. #else
  183. double floor(), ldexp(), polevl(), p1evl();
  184. int isnan(), isfinite();
  185. #endif
  186. extern double MAXNUM;
  187. #ifdef INFINITIES
  188. extern double INFINITY;
  189. #endif
  190. double exp10(x) double x;
  191. {
  192. double px, xx;
  193. short n;
  194. #ifdef NANS
  195. if (isnan(x))
  196. return (x);
  197. #endif
  198. if (x > MAXL10) {
  199. #ifdef INFINITIES
  200. return (INFINITY);
  201. #else
  202. mtherr("exp10", OVERFLOW);
  203. return (MAXNUM);
  204. #endif
  205. }
  206. if (x < -MAXL10) /* Would like to use MINLOG but can't */
  207. {
  208. #ifndef INFINITIES
  209. mtherr("exp10", UNDERFLOW);
  210. #endif
  211. return (0.0);
  212. }
  213. /* Express 10**x = 10**g 2**n
  214. * = 10**g 10**( n log10(2) )
  215. * = 10**( g + n log10(2) )
  216. */
  217. px = floor(LOG210 * x + 0.5);
  218. n = px;
  219. x -= px * LG102A;
  220. x -= px * LG102B;
  221. /* rational approximation for exponential
  222. * of the fractional part:
  223. * 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
  224. */
  225. xx = x * x;
  226. px = x * polevl(xx, P, 3);
  227. x = px / (p1evl(xx, Q, 3) - px);
  228. x = 1.0 + ldexp(x, 1);
  229. /* multiply by power of 2 */
  230. x = ldexp(x, n);
  231. return (x);
  232. }