lib/math/exp.c

   1 /*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7 /* $Header$ */
   8
   9 #include        <math.h>
  10 #include        <float.h>
  11 #include        <errno.h>
  12 #include        "localmath.h"
  13
  14
  15 double
  16 exp(double x)
  17 {
  18         /*      Algorithm and coefficients from:
  19                         "Software manual for the elementary functions"
  20                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  21         */
  22
  23         static double p[] = {
  24                 0.25000000000000000000e+0,
  25                 0.75753180159422776666e-2,
  26                 0.31555192765684646356e-4
  27         };
  28
  29         static double q[] = {
  30                 0.50000000000000000000e+0,
  31                 0.56817302698551221787e-1,
  32                 0.63121894374398503557e-3,
  33                 0.75104028399870046114e-6
  34         };
  35         double  xn, g;
  36         int     n;
  37         int     negative = x < 0;
  38
  39         if (__IsNan(x)) {
  40                 errno = EDOM;
  41                 return x;
  42         }
  43         if (x < M_LN_MIN_D) {
  44                 errno = ERANGE;
  45                 return 0.0;
  46         }
  47         if (x > M_LN_MAX_D) {
  48                 errno = ERANGE;
  49                 return HUGE_VAL;
  50         }
  51
  52         if (negative) x = -x;
  53
  54         /* ??? avoid underflow ??? */
  55
  56         n = x * M_LOG2E + 0.5;  /* 1/ln(2) = log2(e), 0.5 added for rounding */
  57         xn = n;
  58         {
  59                 double  x1 = (long) x;
  60                 double  x2 = x - x1;
  61
  62                 g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
  63         }
  64         if (negative) {
  65                 g = -g;
  66                 n = -n;
  67         }
  68         xn = g * g;
  69         x = g * POLYNOM2(xn, p);
  70         n += 1;
  71         return (ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
  72 }