expl.c 3.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127
  1. /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expl.c */
  2. /*
  3. * Copyright (c) 2008 Stephen L. Moshier <[email protected]>
  4. *
  5. * Permission to use, copy, modify, and distribute this software for any
  6. * purpose with or without fee is hereby granted, provided that the above
  7. * copyright notice and this permission notice appear in all copies.
  8. *
  9. * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  10. * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  11. * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  12. * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  13. * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  14. * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  15. * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  16. */
  17. /*
  18. * Exponential function, long double precision
  19. *
  20. *
  21. * SYNOPSIS:
  22. *
  23. * long double x, y, expl();
  24. *
  25. * y = expl( x );
  26. *
  27. *
  28. * DESCRIPTION:
  29. *
  30. * Returns e (2.71828...) raised to the x power.
  31. *
  32. * Range reduction is accomplished by separating the argument
  33. * into an integer k and fraction f such that
  34. *
  35. * x k f
  36. * e = 2 e.
  37. *
  38. * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
  39. * in the basic range [-0.5 ln 2, 0.5 ln 2].
  40. *
  41. *
  42. * ACCURACY:
  43. *
  44. * Relative error:
  45. * arithmetic domain # trials peak rms
  46. * IEEE +-10000 50000 1.12e-19 2.81e-20
  47. *
  48. *
  49. * Error amplification in the exponential function can be
  50. * a serious matter. The error propagation involves
  51. * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
  52. * which shows that a 1 lsb error in representing X produces
  53. * a relative error of X times 1 lsb in the function.
  54. * While the routine gives an accurate result for arguments
  55. * that are exactly represented by a long double precision
  56. * computer number, the result contains amplified roundoff
  57. * error for large arguments not exactly represented.
  58. *
  59. *
  60. * ERROR MESSAGES:
  61. *
  62. * message condition value returned
  63. * exp underflow x < MINLOG 0.0
  64. * exp overflow x > MAXLOG MAXNUM
  65. *
  66. */
  67. #include "libm.h"
  68. #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
  69. long double expl(long double x)
  70. {
  71. return exp(x);
  72. }
  73. #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
  74. static const long double P[3] = {
  75. 1.2617719307481059087798E-4L,
  76. 3.0299440770744196129956E-2L,
  77. 9.9999999999999999991025E-1L,
  78. };
  79. static const long double Q[4] = {
  80. 3.0019850513866445504159E-6L,
  81. 2.5244834034968410419224E-3L,
  82. 2.2726554820815502876593E-1L,
  83. 2.0000000000000000000897E0L,
  84. };
  85. static const long double
  86. C1 = 6.9314575195312500000000E-1L,
  87. C2 = 1.4286068203094172321215E-6L,
  88. MAXLOGL = 1.1356523406294143949492E4L,
  89. MINLOGL = -1.13994985314888605586758E4L,
  90. LOG2EL = 1.4426950408889634073599E0L;
  91. long double expl(long double x)
  92. {
  93. long double px, xx;
  94. int n;
  95. if (isnan(x))
  96. return x;
  97. if (x > MAXLOGL)
  98. return INFINITY;
  99. if (x < MINLOGL)
  100. return 0.0;
  101. /* Express e**x = e**g 2**n
  102. * = e**g e**(n loge(2))
  103. * = e**(g + n loge(2))
  104. */
  105. px = floorl(LOG2EL * x + 0.5); /* floor() truncates toward -infinity. */
  106. n = px;
  107. x -= px * C1;
  108. x -= px * C2;
  109. /* rational approximation for exponential
  110. * of the fractional part:
  111. * e**x = 1 + 2x P(x**2)/(Q(x**2) - P(x**2))
  112. */
  113. xx = x * x;
  114. px = x * __polevll(xx, P, 2);
  115. x = px/(__polevll(xx, Q, 3) - px);
  116. x = 1.0 + 2.0 * x;
  117. x = scalbnl(x, n);
  118. return x;
  119. }
  120. #endif