home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ieee / eexp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.2 KB  |  60 lines
  1. /*                            xexp.c        */
  2. /* exponential function check routine */
  3. /* by Stephen L. Moshier. */
  4.  
  5.  
  6. #include "ehead.h"
  7.  
  8.  
  9. void eexp( x, y )
  10. unsigned short *x, *y;
  11. {
  12. unsigned short num[NE], den[NE], x2[NE];
  13. long i;
  14. unsigned short sign;
  15.  
  16. /* range reduction theory: x = i + f, 0<=f<1;
  17.  * e**x = e**i * e**f 
  18.  * e**i = 2**(i/log 2).
  19.  * Let i/log2 = i1 + f1, 0<=f1<1.
  20.  * Then e**i = 2**i1 * 2**f1, so
  21.  * e**x = 2**i1 * e**(log 2 * f1) * e**f.
  22.  */
  23.  
  24. if( (x[NE-1] & (unsigned short )0x7fff) == 0 )
  25.     {
  26.     emov( eone, y );
  27.     return;
  28.     }
  29. emov(x, x2);
  30. sign = x2[NE-1] & 0x8000;
  31. if( sign )
  32.     x2[NE-1] &= 0x7fff;
  33. eifrac( x2, &i, num );        /* x = i + f        */
  34.  
  35. if( i != 0 )
  36.  {
  37.  ltoe( &i, den );        /* floating point i    */
  38.  ediv( elog2, den, den );    /* i/log 2        */
  39.  eifrac( den, &i, den );    /* i/log 2  =  i1 + f1    */
  40.  emul( elog2, den, den );    /* log 2 * f1        */
  41.  eadd( den, num, x2 );        /* log 2 * f1  + f    */
  42.  }
  43.  
  44. /*x2[NE-1] -= 1;*/
  45. eldexp( x2, -1L, x2 ); /* divide by 2 */
  46. etanh( x2, x2 );    /* tanh( x/2 )            */
  47. eadd( x2, eone, num );    /* 1 + tanh            */
  48. eneg( x2 );
  49. eadd( x2, eone, den );    /* 1 - tanh            */
  50. ediv( den, num, y );    /* (1 + tanh)/(1 - tanh)    */
  51.  
  52. /*y[NE-1] += i;*/
  53. if( sign )
  54.     {
  55.     ediv( y, eone, y );
  56.     i = -i;
  57.     }
  58. eldexp( y, i, y );    /* multiply by 2**i */
  59. }
  60.