exp2f.c 1.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116
  1. /* exp2f.c
  2. *
  3. * Base 2 exponential function
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * float x, y, exp2f();
  10. *
  11. * y = exp2f( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns 2 raised to the x power.
  18. *
  19. * Range reduction is accomplished by separating the argument
  20. * into an integer k and fraction f such that
  21. * x k f
  22. * 2 = 2 2.
  23. *
  24. * A polynomial approximates 2**x in the basic range [-0.5, 0.5].
  25. *
  26. *
  27. * ACCURACY:
  28. *
  29. * Relative error:
  30. * arithmetic domain # trials peak rms
  31. * IEEE -127,+127 100000 1.7e-7 2.8e-8
  32. *
  33. *
  34. * See exp.c for comments on error amplification.
  35. *
  36. *
  37. * ERROR MESSAGES:
  38. *
  39. * message condition value returned
  40. * exp underflow x < -MAXL2 0.0
  41. * exp overflow x > MAXL2 MAXNUMF
  42. *
  43. * For IEEE arithmetic, MAXL2 = 127.
  44. */
  45. /*
  46. Cephes Math Library Release 2.2: June, 1992
  47. Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
  48. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  49. */
  50. #include <math.h>
  51. static char fname[] = {"exp2f"};
  52. static float P[] = {
  53. 1.535336188319500E-004,
  54. 1.339887440266574E-003,
  55. 9.618437357674640E-003,
  56. 5.550332471162809E-002,
  57. 2.402264791363012E-001,
  58. 6.931472028550421E-001
  59. };
  60. #define MAXL2 127.0
  61. #define MINL2 -127.0
  62. extern float MAXNUMF;
  63. float polevlf(float, float *, int), floorf(float), ldexpf(float, int);
  64. float exp2f( float xx )
  65. {
  66. float x, px;
  67. int i0;
  68. x = xx;
  69. if( x > MAXL2)
  70. {
  71. mtherr( fname, OVERFLOW );
  72. return( MAXNUMF );
  73. }
  74. if( x < MINL2 )
  75. {
  76. mtherr( fname, UNDERFLOW );
  77. return(0.0);
  78. }
  79. /* The following is necessary because range reduction blows up: */
  80. if( x == 0 )
  81. return(1.0);
  82. /* separate into integer and fractional parts */
  83. px = floorf(x);
  84. i0 = px;
  85. x = x - px;
  86. if( x > 0.5 )
  87. {
  88. i0 += 1;
  89. x -= 1.0;
  90. }
  91. /* rational approximation
  92. * exp2(x) = 1.0 + xP(x)
  93. */
  94. px = 1.0 + x * polevlf( x, P, 5 );
  95. /* scale by power of 2 */
  96. px = ldexpf( px, i0 );
  97. return(px);
  98. }