libm/float/exp10f.c

   1 /*                                                      exp10f.c
   2  *
   3  *      Base 10 exponential function
   4  *      (Common antilogarithm)
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * float x, y, exp10f();
  11  *
  12  * y = exp10f( x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns 10 raised to the x power.
  19  *
  20  * Range reduction is accomplished by expressing the argument
  21  * as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
  22  * A polynomial approximates 10**f.
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *                      Relative error:
  29  * arithmetic   domain     # trials      peak         rms
  30  *    IEEE      -38,+38     100000      9.8e-8      2.8e-8
  31  *
  32  * ERROR MESSAGES:
  33  *
  34  *   message         condition      value returned
  35  * exp10 underflow    x < -MAXL10        0.0
  36  * exp10 overflow     x > MAXL10       MAXNUM
  37  *
  38  * IEEE single arithmetic: MAXL10 = 38.230809449325611792.
  39  *
  40  */
  41 \f
  42 /*
  43 Cephes Math Library Release 2.2:  June, 1992
  44 Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
  45 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  46 */
  47
  48
  49 #include <math.h>
  50
  51 static float P[] = {
  52  2.063216740311022E-001,
  53  5.420251702225484E-001,
  54  1.171292686296281E+000,
  55  2.034649854009453E+000,
  56  2.650948748208892E+000,
  57  2.302585167056758E+000
  58 };
  59
  60 /*static float LOG102 = 3.01029995663981195214e-1;*/
  61 static float LOG210 = 3.32192809488736234787e0;
  62 static float LG102A = 3.00781250000000000000E-1;
  63 static float LG102B = 2.48745663981195213739E-4;
  64 static float MAXL10 = 38.230809449325611792;
  65
  66
  67
  68
  69 extern float MAXNUMF;
  70
  71 float floorf(float), ldexpf(float, int), polevlf(float, float *, int);
  72
  73 float exp10f(float xx)
  74 {
  75 float x, px, qx;
  76 short n;
  77
  78 x = xx;
  79 if( x > MAXL10 )
  80         {
  81         mtherr( "exp10f", OVERFLOW );
  82         return( MAXNUMF );
  83         }
  84
  85 if( x < -MAXL10 )       /* Would like to use MINLOG but can't */
  86         {
  87         mtherr( "exp10f", UNDERFLOW );
  88         return(0.0);
  89         }
  90
  91 /* The following is necessary because range reduction blows up: */
  92 if( x == 0 )
  93         return(1.0);
  94
  95 /* Express 10**x = 10**g 2**n
  96  *   = 10**g 10**( n log10(2) )
  97  *   = 10**( g + n log10(2) )
  98  */
  99 px = x * LOG210;
 100 qx = floorf( px + 0.5 );
 101 n = qx;
 102 x -= qx * LG102A;
 103 x -= qx * LG102B;
 104
 105 /* rational approximation for exponential
 106  * of the fractional part:
 107  * 10**x - 1  =  2x P(x**2)/( Q(x**2) - P(x**2) )
 108  */
 109 px = 1.0 + x * polevlf( x, P, 5 );
 110
 111 /* multiply by power of 2 */
 112 x = ldexpf( px, n );
 113
 114 return(x);
 115 }