a) you can do a finite sum evaluation (with n
exp(x) = sum( n=0..infinity, (x^n)/(n!) );

  double exp_approx( double x, unsigned int n ){  double result = 0;  double fact = 1;  while( n > 0 )  {    result = (fact + result) * x;    fact *= n--;  }  result /= fact;  return 1.0+result;}int main(int argc, char* argv[]){  int i;  double x = 0;  for( i = 0; i < 100; ++ i)  {    double e1, e2, error = 0;    while( error <= 1e-6 )    {      e1 = exp_approx( x, i );      e2 = exp( x );      printf( "%6.3f %.20lf\r", x, error = fabs(e1-e2)/e2 );      x += 0.001;    }    printf( "\n" );  }  return 0;}  

Or, alternatively, try a decomposition of the IEEE float and a lookup table.

(not quite tested code)

  double exp_lookup_table[256];void exp_lookup_init(){   int i;   for( i = 0; i < 256; ++i )     exp_lookup_table[i] = expf( 1.0 + i / 256.0 ); // won''t work with denormalized numbers;}float exp_lookup( float in ) // yes, we take and return float not double.{  typedef union { int i; float f; } f2i_t;  f2i_t num;  num.f = in  bool negative = (num & 0x80000000) != 0;  long exponent = ((num & 0x7F800000) >> 23) - 127;  long mantissa = ((num & 0x007F8000) >> 15); // we''ll use a 256 entry lookup table  double result = exp_lookup_table[mantissa]; // yep it''s a double.  // Now, unfortunately, we have processing to do...  while( exponent > 0 )   {     result *= result;    ++exponent;  }  while( exponent < 0 )   {     result = sqrtf( result ); // THIS WILL BE SLOW    ++exponent;  }    if( negative ) result = 1.0 / result;    return (float)result;}  

Most of the ideas you will find on the www are based on a fast 2 power n computation and limited developpement of exp(x)

exp(x) = 1+x+x2/2+x3/6 + xn/n! + ...

But I have a quicker solution.

Computing 2^n is simply setting the exponent of a floating point number. Read about the IEEE standard floating point representation.

Then as exp(x) = 2^(x/ln(2)), compute y=x*(1/ln(2))

Lets define as i=int(y), f=y-int(y)

Thus exp(x) = 2^i * 2^y. y beeing between 0 and 1 you can use the fast fixed point conversion trick (adding a big number) and it will be an index into your 2^y (0
// 16 bit precision exponential
// This works on any PC. It should work on most machines with some small precautions. (#defines)

float exp_16b( float x )
{
byte y[4];
int i,f,g,h; // 32 bits integers

((float*)y)[0] = x*inv_ln2 + magic; // ( store as a 8.16 fixed point )

// Big endian/little endian pb here ? Verify please !
i = y[2]; // integer part.
f = y[1]; // fractionnal part (MSB)
g = y[0]; // fractionnal part (LSB)
h = ((int*)lutExp)[f][0];
h+= i<<23; // Add i to the exponent <=> multiply by 2^i
((int*)y)[0]=h; // store
y*= lutExp[g][1];

return (y);
}

I guess it takes around 10 clock cycles.

You could extend it easilly to 24 or 32 bits by storing a 8.24 or 8.32 fixed point into a double (beware 64 bit alignement then). And 24 bits is the precision of floats !

I dont think anyone can provide a faster code for 16 bits of precision.

Oops : Fruny posted meanwhile. Roughly the same thing but mine is better (add to exponent). Just forgot to check the sign. If x can be negative then exp(-x) = exp(x), thus just add a test of the signbit (y[0]&(1<<31)) and h-=i<<23.


Edited by - Charles B on January 20, 2002 7:00:26 PM

No comparison . Told you ... 10 cycles.