exp.c 1.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172
  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. #include <math.h>
  9. #include <float.h>
  10. #include <errno.h>
  11. #include "localmath.h"
  12. double
  13. exp(double x)
  14. {
  15. /* Algorithm and coefficients from:
  16. "Software manual for the elementary functions"
  17. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  18. */
  19. static double p[] = {
  20. 0.25000000000000000000e+0,
  21. 0.75753180159422776666e-2,
  22. 0.31555192765684646356e-4
  23. };
  24. static double q[] = {
  25. 0.50000000000000000000e+0,
  26. 0.56817302698551221787e-1,
  27. 0.63121894374398503557e-3,
  28. 0.75104028399870046114e-6
  29. };
  30. double xn, g;
  31. int n;
  32. int negative = x < 0;
  33. if (__IsNan(x)) {
  34. errno = EDOM;
  35. return x;
  36. }
  37. if (x < M_LN_MIN_D) {
  38. errno = ERANGE;
  39. return 0.0;
  40. }
  41. if (x > M_LN_MAX_D) {
  42. errno = ERANGE;
  43. return HUGE_VAL;
  44. }
  45. if (negative) x = -x;
  46. /* ??? avoid underflow ??? */
  47. n = x * M_LOG2E + 0.5; /* 1/ln(2) = log2(e), 0.5 added for rounding */
  48. xn = n;
  49. {
  50. double x1 = (long) x;
  51. double x2 = x - x1;
  52. g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
  53. }
  54. if (negative) {
  55. g = -g;
  56. n = -n;
  57. }
  58. xn = g * g;
  59. x = g * POLYNOM2(xn, p);
  60. n += 1;
  61. return (ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
  62. }