lib/libm/exp10.c

   1 /*                                                      exp10.c
   2  *
   3  *      Base 10 exponential function
   4  *      (Common antilogarithm)
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * double x, y, exp10();
  11  *
  12  * y = exp10( 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  * The Pade' form
  23  *
  24  *    1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
  25  *
  26  * is used to approximate 10**f.
  27  *
  28  *
  29  *
  30  * ACCURACY:
  31  *
  32  *                      Relative error:
  33  * arithmetic   domain     # trials      peak         rms
  34  *    IEEE     -307,+307    30000       2.2e-16     5.5e-17
  35  * Test result from an earlier version (2.1):
  36  *    DEC       -38,+38     70000       3.1e-17     7.0e-18
  37  *
  38  * ERROR MESSAGES:
  39  *
  40  *   message         condition      value returned
  41  * exp10 underflow    x < -MAXL10        0.0
  42  * exp10 overflow     x > MAXL10       MAXNUM
  43  *
  44  * DEC arithmetic: MAXL10 = 38.230809449325611792.
  45  * IEEE arithmetic: MAXL10 = 308.2547155599167.
  46  *
  47  */
  48 \f
  49 /*
  50 Cephes Math Library Release 2.8:  June, 2000
  51 Copyright 1984, 1991, 2000 by Stephen L. Moshier
  52 */
  53
  54
  55 #include "mconf.h"
  56
  57 #ifdef UNK
  58 static double P[] = {
  59  4.09962519798587023075E-2,
  60  1.17452732554344059015E1,
  61  4.06717289936872725516E2,
  62  2.39423741207388267439E3,
  63 };
  64 static double Q[] = {
  65 /* 1.00000000000000000000E0,*/
  66  8.50936160849306532625E1,
  67  1.27209271178345121210E3,
  68  2.07960819286001865907E3,
  69 };
  70 /* static double LOG102 = 3.01029995663981195214e-1; */
  71 static double LOG210 = 3.32192809488736234787e0;
  72 static double LG102A = 3.01025390625000000000E-1;
  73 static double LG102B = 4.60503898119521373889E-6;
  74 /* static double MAXL10 = 38.230809449325611792; */
  75 static double MAXL10 = 308.2547155599167;
  76 #endif
  77
  78 #ifdef DEC
  79 static unsigned short P[] = {
  80 0037047,0165657,0114061,0067234,
  81 0041073,0166243,0123052,0144643,
  82 0042313,0055720,0024032,0047443,
  83 0043025,0121714,0070232,0050007,
  84 };
  85 static unsigned short Q[] = {
  86 /*0040200,0000000,0000000,0000000,*/
  87 0041652,0027756,0071216,0050075,
  88 0042637,0001367,0077263,0136017,
  89 0043001,0174673,0024157,0133416,
  90 };
  91 /*
  92 static unsigned short L102[] = {0037632,0020232,0102373,0147770};
  93 #define LOG102 *(double *)L102
  94 */
  95 static unsigned short L210[] = {0040524,0115170,0045715,0015613};
  96 #define LOG210 *(double *)L210
  97 static unsigned short L102A[] = {0037632,0020000,0000000,0000000,};
  98 #define LG102A *(double *)L102A
  99 static unsigned short L102B[] = {0033632,0102373,0147767,0114220,};
 100 #define LG102B *(double *)L102B
 101 static unsigned short MXL[] = {0041430,0166131,0047761,0154130,};
 102 #define MAXL10 ( *(double *)MXL )
 103 #endif
 104
 105 #ifdef IBMPC
 106 static unsigned short P[] = {
 107 0x2dd4,0xf306,0xfd75,0x3fa4,
 108 0x5934,0x74c5,0x7d94,0x4027,
 109 0x49e4,0x0503,0x6b7a,0x4079,
 110 0x4a01,0x8e13,0xb479,0x40a2,
 111 };
 112 static unsigned short Q[] = {
 113 /*0x0000,0x0000,0x0000,0x3ff0,*/
 114 0xca08,0xce51,0x45fd,0x4055,
 115 0x7782,0xefd6,0xe05e,0x4093,
 116 0xf6e2,0x650d,0x3f37,0x40a0,
 117 };
 118 /*
 119 static unsigned short L102[] = {0x79ff,0x509f,0x4413,0x3fd3};
 120 #define LOG102 *(double *)L102
 121 */
 122 static unsigned short L210[] = {0xa371,0x0979,0x934f,0x400a};
 123 #define LOG210 *(double *)L210
 124 static unsigned short L102A[] = {0x0000,0x0000,0x4400,0x3fd3,};
 125 #define LG102A *(double *)L102A
 126 static unsigned short L102B[] = {0xf312,0x79fe,0x509f,0x3ed3,};
 127 #define LG102B *(double *)L102B
 128 static double MAXL10 = 308.2547155599167;
 129 #endif
 130
 131 #ifdef MIEEE
 132 static unsigned short P[] = {
 133 0x3fa4,0xfd75,0xf306,0x2dd4,
 134 0x4027,0x7d94,0x74c5,0x5934,
 135 0x4079,0x6b7a,0x0503,0x49e4,
 136 0x40a2,0xb479,0x8e13,0x4a01,
 137 };
 138 static unsigned short Q[] = {
 139 /*0x3ff0,0x0000,0x0000,0x0000,*/
 140 0x4055,0x45fd,0xce51,0xca08,
 141 0x4093,0xe05e,0xefd6,0x7782,
 142 0x40a0,0x3f37,0x650d,0xf6e2,
 143 };
 144 /*
 145 static unsigned short L102[] = {0x3fd3,0x4413,0x509f,0x79ff};
 146 #define LOG102 *(double *)L102
 147 */
 148 static unsigned short L210[] = {0x400a,0x934f,0x0979,0xa371};
 149 #define LOG210 *(double *)L210
 150 static unsigned short L102A[] = {0x3fd3,0x4400,0x0000,0x0000,};
 151 #define LG102A *(double *)L102A
 152 static unsigned short L102B[] = {0x3ed3,0x509f,0x79fe,0xf312,};
 153 #define LG102B *(double *)L102B
 154 static double MAXL10 = 308.2547155599167;
 155 #endif
 156
 157 #ifdef ANSIPROT
 158 extern double floor ( double );
 159 extern double ldexp ( double, int );
 160 extern double polevl ( double, void *, int );
 161 extern double p1evl ( double, void *, int );
 162 extern int isnan ( double );
 163 extern int isfinite ( double );
 164 #else
 165 double floor(), ldexp(), polevl(), p1evl();
 166 int isnan(), isfinite();
 167 #endif
 168 extern double MAXNUM;
 169 #ifdef INFINITIES
 170 extern double INFINITY;
 171 #endif
 172
 173 double exp10(x)
 174 double x;
 175 {
 176 double px, xx;
 177 short n;
 178
 179 #ifdef NANS
 180 if( isnan(x) )
 181         return(x);
 182 #endif
 183 if( x > MAXL10 )
 184         {
 185 #ifdef INFINITIES
 186         return( INFINITY );
 187 #else
 188         mtherr( "exp10", OVERFLOW );
 189         return( MAXNUM );
 190 #endif
 191         }
 192
 193 if( x < -MAXL10 )       /* Would like to use MINLOG but can't */
 194         {
 195 #ifndef INFINITIES
 196         mtherr( "exp10", UNDERFLOW );
 197 #endif
 198         return(0.0);
 199         }
 200
 201 /* Express 10**x = 10**g 2**n
 202  *   = 10**g 10**( n log10(2) )
 203  *   = 10**( g + n log10(2) )
 204  */
 205 px = floor( LOG210 * x + 0.5 );
 206 n = px;
 207 x -= px * LG102A;
 208 x -= px * LG102B;
 209
 210 /* rational approximation for exponential
 211  * of the fractional part:
 212  * 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
 213  */
 214 xx = x * x;
 215 px = x * polevl( xx, P, 3 );
 216 x =  px/( p1evl( xx, Q, 3 ) - px );
 217 x = 1.0 + ldexp( x, 1 );
 218
 219 /* multiply by power of 2 */
 220 x = ldexp( x, n );
 221
 222 return(x);
 223 }