sys/src/libc/port/exp.c

   1 /*
   2         exp returns the exponential function of its
   3         floating-point argument.
   4
   5         The coefficients are #1069 from Hart and Cheney. (22.35D)
   6 */
   7
   8 #include <u.h>
   9 #include <libc.h>
  10
  11 #define p0  .2080384346694663001443843411e7
  12 #define p1  .3028697169744036299076048876e5
  13 #define p2  .6061485330061080841615584556e2
  14 #define q0  .6002720360238832528230907598e7
  15 #define q1  .3277251518082914423057964422e6
  16 #define q2  .1749287689093076403844945335e4
  17 #define log2e  1.4426950408889634073599247
  18 #define sqrt2  1.4142135623730950488016887
  19 #define maxf  10000
  20
  21 double
  22 exp(double arg)
  23 {
  24         double fract, temp1, temp2, xsq;
  25         int ent;
  26
  27         if(arg == 0)
  28                 return 1;
  29         if(arg < -maxf)
  30                 return 0;
  31         if(arg > maxf)
  32                 return Inf(1);
  33         arg *= log2e;
  34         ent = floor(arg);
  35         fract = (arg-ent) - 0.5;
  36         xsq = fract*fract;
  37         temp1 = ((p2*xsq+p1)*xsq+p0)*fract;
  38         temp2 = ((xsq+q2)*xsq+q1)*xsq + q0;
  39         return ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent);
  40 }