lib/ack/libp/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
   8 /* $Header$ */
   9 #define __NO_DEFS
  10 #include <math.h>
  11 #include <pc_err.h>
  12 extern  _trp();
  13
  14 #if __STDC__
  15 #include <float.h>
  16 #include <pc_math.h>
  17 #define M_MIN_D DBL_MIN
  18 #define M_MAX_D DBL_MAX
  19 #define M_DMINEXP DBL_MIN_EXP
  20 #endif
  21 #undef HUGE
  22 #define HUGE    1e1000
  23
  24 static double
  25 Ldexp(fl,exp)
  26         double fl;
  27         int exp;
  28 {
  29         extern double _fef();
  30         int sign = 1;
  31         int currexp;
  32
  33         if (fl<0) {
  34                 fl = -fl;
  35                 sign = -1;
  36         }
  37         fl = _fef(fl,&currexp);
  38         exp += currexp;
  39         if (exp > 0) {
  40                 while (exp>30) {
  41                         fl *= (double) (1L << 30);
  42                         exp -= 30;
  43                 }
  44                 fl *= (double) (1L << exp);
  45         }
  46         else    {
  47                 while (exp<-30) {
  48                         fl /= (double) (1L << 30);
  49                         exp += 30;
  50                 }
  51                 fl /= (double) (1L << -exp);
  52         }
  53         return sign * fl;
  54 }
  55
  56 double
  57 _exp(x)
  58         double x;
  59 {
  60         /*      Algorithm and coefficients from:
  61                         "Software manual for the elementary functions"
  62                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  63         */
  64
  65         static double p[] = {
  66                 0.25000000000000000000e+0,
  67                 0.75753180159422776666e-2,
  68                 0.31555192765684646356e-4
  69         };
  70
  71         static double q[] = {
  72                 0.50000000000000000000e+0,
  73                 0.56817302698551221787e-1,
  74                 0.63121894374398503557e-3,
  75                 0.75104028399870046114e-6
  76         };
  77         double  xn, g;
  78         int     n;
  79         int     negative = x < 0;
  80
  81         if (x <= M_LN_MIN_D) {
  82                 g = M_MIN_D/4.0;
  83
  84                 if (g != 0.0) {
  85                         /* unnormalized numbers apparently exist */
  86                         if (x < (M_LN2 * (M_DMINEXP - 53))) return 0.0;
  87                 }
  88                 else {
  89                         if (x < M_LN_MIN_D) return 0.0;
  90                         return M_MIN_D;
  91                 }
  92         }
  93         if (x >= M_LN_MAX_D) {
  94                 if (x > M_LN_MAX_D) {
  95                         _trp(EEXP);
  96                         return HUGE;
  97                 }
  98                 return M_MAX_D;
  99         }
 100         if (negative) x = -x;
 101
 102         n = x * M_LOG2E + 0.5;  /* 1/ln(2) = log2(e), 0.5 added for rounding */
 103         xn = n;
 104         {
 105                 double  x1 = (long) x;
 106                 double  x2 = x - x1;
 107
 108                 g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
 109         }
 110         if (negative) {
 111                 g = -g;
 112                 n = -n;
 113         }
 114         xn = g * g;
 115         x = g * POLYNOM2(xn, p);
 116         n += 1;
 117         return (Ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
 118 }