expx2.c 1.5 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091
  1. /* expx2.c
  2. *
  3. * Exponential of squared argument
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, y, expx2();
  10. * int sign;
  11. *
  12. * y = expx2( x, sign );
  13. *
  14. *
  15. *
  16. * DESCRIPTION:
  17. *
  18. * Computes y = exp(x*x) while suppressing error amplification
  19. * that would ordinarily arise from the inexactness of the
  20. * exponential argument x*x.
  21. *
  22. * If sign < 0, the result is inverted; i.e., y = exp(-x*x) .
  23. *
  24. *
  25. * ACCURACY:
  26. *
  27. * Relative error:
  28. * arithmetic domain # trials peak rms
  29. * IEEE -26.6, 26.6 10^7 3.9e-16 8.9e-17
  30. *
  31. */
  32. /*
  33. Cephes Math Library Release 2.9: June, 2000
  34. Copyright 2000 by Stephen L. Moshier
  35. */
  36. #include "mconf.h"
  37. #ifdef ANSIPROT
  38. extern double fabs(double);
  39. extern double floor(double);
  40. extern double exp(double);
  41. #else
  42. double fabs();
  43. double floor();
  44. double exp();
  45. #endif
  46. #ifdef DEC
  47. #define M 32.0
  48. #define MINV .03125
  49. #else
  50. #define M 128.0
  51. #define MINV .0078125
  52. #endif
  53. extern double MAXLOG;
  54. extern double INFINITY;
  55. double expx2(x, sign) double x;
  56. int sign;
  57. {
  58. double u, u1, m, f;
  59. x = fabs(x);
  60. if (sign < 0)
  61. x = -x;
  62. /* Represent x as an exact multiple of M plus a residual.
  63. M is a power of 2 chosen so that exp(m * m) does not overflow
  64. or underflow and so that |x - m| is small. */
  65. m = MINV * floor(M * x + 0.5);
  66. f = x - m;
  67. /* x^2 = m^2 + 2mf + f^2 */
  68. u = m * m;
  69. u1 = 2 * m * f + f * f;
  70. if (sign < 0) {
  71. u = -u;
  72. u1 = -u1;
  73. }
  74. if ((u + u1) > MAXLOG)
  75. return (INFINITY);
  76. /* u is exact, u1 is small. */
  77. u = exp(u) * exp(u1);
  78. return (u);
  79. }