home *** CD-ROM | disk | FTP | other *** search
/ Tricks of the Windows Gam…ming Gurus (2nd Edition) / Disc2.iso / vc98 / crt / src / xlexp.c < prev    next >
C/C++ Source or Header  |  1998-06-16  |  2KB  |  65 lines
  1. /* _LExp function */
  2. #include "wctype.h"
  3. #include "xmath.h"
  4. _STD_BEGIN
  5.  
  6. /* coefficients */
  7. static const long double p[] = {    /* courtesy Dr. Tim Prince */
  8.     42.038913947607355L,
  9.     10096.353102778762831L,
  10.     333228.767219512631062L};
  11. static const long double q[] = {    /* courtesy Dr. Tim Prince */
  12.     1.0L,
  13.     841.167880526530790L,
  14.     75730.834075476293976L,
  15.     666457.534439025262146L};
  16. static const long double c1 = 22713.0 / 32768.0;
  17. static const long double c2 = 1.428606820309417232e-6L;
  18. static const long double hugexp = LHUGE_EXP;
  19. static const long double invln2 = 1.4426950408889634074L;
  20.  
  21. _CRTIMP2 short __cdecl _LExp(long double *px, long double y, short eoff)
  22.     {    /* compute y*e^(*px), (*px) finite, |y| not huge */
  23.     if (*px < -hugexp || y == 0)
  24.         {    /* certain underflow */
  25.         *px = 0;
  26.         return (0);
  27.         }
  28.     else if (hugexp < *px)
  29.         {    /* certain overflow */
  30.         *px = _LInf._L;
  31.         return (INF);
  32.         }
  33.     else
  34.         {    /* xexp won't overflow */
  35.         long double g = *px * invln2;
  36.         short xexp = g + (g < 0 ? - 0.5 : + 0.5);
  37.  
  38.         g = xexp;
  39.         g = (*px - g * c1) - g * c2;
  40.         if (-_LEps._L < g && g < _LEps._L)
  41.             *px = y;
  42.         else
  43.             {    /* g*g worth computing */
  44.             const long double z = g * g;
  45.             const long double w = ((z + q[1]) * z + q[2]) * z
  46.                 + q[3];
  47.  
  48.             g *= (p[0] * z + p[1]) * z + p[2];
  49.             *px = (w + g) / (w - g) * 2 * y;
  50.             --xexp;
  51.             }
  52.         return (_LDscale(px, (long)xexp + eoff));
  53.         }
  54.     }
  55. _STD_END
  56.  
  57. /*
  58.  * Copyright (c) 1994 by P.J. Plauger.  ALL RIGHTS RESERVED. 
  59.  * Consult your license regarding permissions and restrictions.
  60.  */
  61.  
  62. /*
  63. 941029 pjp: added _STD machinery
  64.  */
  65.