e_expl.c 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145
  1. /* $OpenBSD: e_expl.c,v 1.3 2013/11/12 20:35:18 martynas Exp $ */
  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. /* expl.c
  18. *
  19. * Exponential function, 128-bit long double precision
  20. *
  21. *
  22. *
  23. * SYNOPSIS:
  24. *
  25. * long double x, y, expl();
  26. *
  27. * y = expl( x );
  28. *
  29. *
  30. *
  31. * DESCRIPTION:
  32. *
  33. * Returns e (2.71828...) raised to the x power.
  34. *
  35. * Range reduction is accomplished by separating the argument
  36. * into an integer k and fraction f such that
  37. *
  38. * x k f
  39. * e = 2 e.
  40. *
  41. * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
  42. * in the basic range [-0.5 ln 2, 0.5 ln 2].
  43. *
  44. *
  45. * ACCURACY:
  46. *
  47. * Relative error:
  48. * arithmetic domain # trials peak rms
  49. * IEEE +-MAXLOG 100,000 2.6e-34 8.6e-35
  50. *
  51. *
  52. * Error amplification in the exponential function can be
  53. * a serious matter. The error propagation involves
  54. * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
  55. * which shows that a 1 lsb error in representing X produces
  56. * a relative error of X times 1 lsb in the function.
  57. * While the routine gives an accurate result for arguments
  58. * that are exactly represented by a long double precision
  59. * computer number, the result contains amplified roundoff
  60. * error for large arguments not exactly represented.
  61. *
  62. *
  63. * ERROR MESSAGES:
  64. *
  65. * message condition value returned
  66. * exp underflow x < MINLOG 0.0
  67. * exp overflow x > MAXLOG MAXNUM
  68. *
  69. */
  70. /* Exponential function */
  71. #include <float.h>
  72. #include <openlibm_math.h>
  73. #include "math_private.h"
  74. /* Pade' coefficients for exp(x) - 1
  75. Theoretical peak relative error = 2.2e-37,
  76. relative peak error spread = 9.2e-38
  77. */
  78. static long double P[5] = {
  79. 3.279723985560247033712687707263393506266E-10L,
  80. 6.141506007208645008909088812338454698548E-7L,
  81. 2.708775201978218837374512615596512792224E-4L,
  82. 3.508710990737834361215404761139478627390E-2L,
  83. 9.999999999999999999999999999999999998502E-1L
  84. };
  85. static long double Q[6] = {
  86. 2.980756652081995192255342779918052538681E-12L,
  87. 1.771372078166251484503904874657985291164E-8L,
  88. 1.504792651814944826817779302637284053660E-5L,
  89. 3.611828913847589925056132680618007270344E-3L,
  90. 2.368408864814233538909747618894558968880E-1L,
  91. 2.000000000000000000000000000000000000150E0L
  92. };
  93. /* C1 + C2 = ln 2 */
  94. static const long double C1 = -6.93145751953125E-1L;
  95. static const long double C2 = -1.428606820309417232121458176568075500134E-6L;
  96. static const long double LOG2EL = 1.442695040888963407359924681001892137426646L;
  97. static const long double MAXLOGL = 1.1356523406294143949491931077970764891253E4L;
  98. static const long double MINLOGL = -1.143276959615573793352782661133116431383730e4L;
  99. static const long double huge = 0x1p10000L;
  100. #if 0 /* XXX Prevent gcc from erroneously constant folding this. */
  101. static const long double twom10000 = 0x1p-10000L;
  102. #else
  103. static volatile long double twom10000 = 0x1p-10000L;
  104. #endif
  105. long double
  106. expl(long double x)
  107. {
  108. long double px, xx;
  109. int n;
  110. if( x > MAXLOGL)
  111. return (huge*huge); /* overflow */
  112. if( x < MINLOGL )
  113. return (twom10000*twom10000); /* underflow */
  114. /* Express e**x = e**g 2**n
  115. * = e**g e**( n loge(2) )
  116. * = e**( g + n loge(2) )
  117. */
  118. px = floorl( LOG2EL * x + 0.5L ); /* floor() truncates toward -infinity. */
  119. n = px;
  120. x += px * C1;
  121. x += px * C2;
  122. /* rational approximation for exponential
  123. * of the fractional part:
  124. * e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
  125. */
  126. xx = x * x;
  127. px = x * __polevll( xx, P, 4 );
  128. xx = __polevll( xx, Q, 5 );
  129. x = px/( xx - px );
  130. x = 1.0L + x + x;
  131. x = ldexpl( x, n );
  132. return(x);
  133. }