exp.c 1.7 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697
  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. static double
  11. ldexp(fl,exp)
  12. double fl;
  13. int exp;
  14. {
  15. extern double _fef();
  16. int sign = 1;
  17. int currexp;
  18. if (fl<0) {
  19. fl = -fl;
  20. sign = -1;
  21. }
  22. fl = _fef(fl,&currexp);
  23. exp += currexp;
  24. if (exp > 0) {
  25. while (exp>30) {
  26. fl *= (double) (1L << 30);
  27. exp -= 30;
  28. }
  29. fl *= (double) (1L << exp);
  30. }
  31. else {
  32. while (exp<-30) {
  33. fl /= (double) (1L << 30);
  34. exp += 30;
  35. }
  36. fl /= (double) (1L << -exp);
  37. }
  38. return sign * fl;
  39. }
  40. double
  41. _exp(x)
  42. double x;
  43. {
  44. /* Algorithm and coefficients from:
  45. "Software manual for the elementary functions"
  46. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  47. */
  48. static double p[] = {
  49. 0.25000000000000000000e+0,
  50. 0.75753180159422776666e-2,
  51. 0.31555192765684646356e-4
  52. };
  53. static double q[] = {
  54. 0.50000000000000000000e+0,
  55. 0.56817302698551221787e-1,
  56. 0.63121894374398503557e-3,
  57. 0.75104028399870046114e-6
  58. };
  59. double xn, g;
  60. int n;
  61. int negative = x < 0;
  62. if (x <= M_LN_MIN_D) {
  63. return M_MIN_D;
  64. }
  65. if (x >= M_LN_MAX_D) {
  66. if (x > M_LN_MAX_D) error(3);
  67. return M_MAX_D;
  68. }
  69. if (negative) x = -x;
  70. /* ??? avoid underflow ??? */
  71. n = x * M_LOG2E + 0.5; /* 1/ln(2) = log2(e), 0.5 added for rounding */
  72. xn = n;
  73. {
  74. double x1 = (long) x;
  75. double x2 = x - x1;
  76. g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
  77. }
  78. if (negative) {
  79. g = -g;
  80. n = -n;
  81. }
  82. xn = g * g;
  83. x = g * POLYNOM2(xn, p);
  84. n += 1;
  85. return (ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
  86. }