exp.c 2.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118
  1. /*
  2. * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  3. * See the copyright notice in the ACK home directory, in the file "Copyright".
  4. *
  5. * Author: Ceriel J.H. Jacobs
  6. */
  7. /* $Id$ */
  8. #define __NO_DEFS
  9. #include <math.h>
  10. #include <pc_err.h>
  11. extern _trp();
  12. #if __STDC__
  13. #include <float.h>
  14. #include <pc_math.h>
  15. #define M_MIN_D DBL_MIN
  16. #define M_MAX_D DBL_MAX
  17. #define M_DMINEXP DBL_MIN_EXP
  18. #endif
  19. #undef HUGE
  20. #define HUGE 1e1000
  21. static double
  22. Ldexp(fl,exp)
  23. double fl;
  24. int exp;
  25. {
  26. extern double _fef();
  27. int sign = 1;
  28. int currexp;
  29. if (fl<0) {
  30. fl = -fl;
  31. sign = -1;
  32. }
  33. fl = _fef(fl,&currexp);
  34. exp += currexp;
  35. if (exp > 0) {
  36. while (exp>30) {
  37. fl *= (double) (1L << 30);
  38. exp -= 30;
  39. }
  40. fl *= (double) (1L << exp);
  41. }
  42. else {
  43. while (exp<-30) {
  44. fl /= (double) (1L << 30);
  45. exp += 30;
  46. }
  47. fl /= (double) (1L << -exp);
  48. }
  49. return sign * fl;
  50. }
  51. double
  52. _exp(x)
  53. double x;
  54. {
  55. /* Algorithm and coefficients from:
  56. "Software manual for the elementary functions"
  57. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  58. */
  59. static double p[] = {
  60. 0.25000000000000000000e+0,
  61. 0.75753180159422776666e-2,
  62. 0.31555192765684646356e-4
  63. };
  64. static double q[] = {
  65. 0.50000000000000000000e+0,
  66. 0.56817302698551221787e-1,
  67. 0.63121894374398503557e-3,
  68. 0.75104028399870046114e-6
  69. };
  70. double xn, g;
  71. int n;
  72. int negative = x < 0;
  73. if (x <= M_LN_MIN_D) {
  74. g = M_MIN_D/4.0;
  75. if (g != 0.0) {
  76. /* unnormalized numbers apparently exist */
  77. if (x < (M_LN2 * (M_DMINEXP - 53))) return 0.0;
  78. }
  79. else {
  80. if (x < M_LN_MIN_D) return 0.0;
  81. return M_MIN_D;
  82. }
  83. }
  84. if (x >= M_LN_MAX_D) {
  85. if (x > M_LN_MAX_D) {
  86. _trp(EEXP);
  87. return HUGE;
  88. }
  89. return M_MAX_D;
  90. }
  91. if (negative) x = -x;
  92. n = x * M_LOG2E + 0.5; /* 1/ln(2) = log2(e), 0.5 added for rounding */
  93. xn = n;
  94. {
  95. double x1 = (long) x;
  96. double x2 = x - x1;
  97. g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
  98. }
  99. if (negative) {
  100. g = -g;
  101. n = -n;
  102. }
  103. xn = g * g;
  104. x = g * POLYNOM2(xn, p);
  105. n += 1;
  106. return (Ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
  107. }