Incorporated from BSD 4.4-Lite.

diff --git a/sysdeps/generic/exp.c b/sysdeps/generic/exp.c
index dc1b944..9b4f045 100644 (file)
--- a/sysdeps/generic/exp.c
+++ b/sysdeps/generic/exp.c
@@ -1,6 +1,6 @@
 /*
- * Copyright (c) 1985 Regents of the University of California.
- * All rights reserved.
+ * Copyright (c) 1985, 1993
+ *     The Regents of the University of California.  All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions
@@ -32,7 +32,7 @@
  */
 
 #ifndef lint
-static char sccsid[] = "@(#)exp.c      5.6 (Berkeley) 10/9/90";
+static char sccsid[] = "@(#)exp.c      8.1 (Berkeley) 6/4/93";
 #endif /* not lint */
 
 /* EXP(X)
@@ -156,3 +156,48 @@ double x;
        /* exp(INF) is INF, exp(+big#) overflows to INF */
            return( finite(x) ?  scalb(1.0,5000)  : x);
 }
+
+/* returns exp(r = x + c) for |c| < |x| with no overlap.  */
+
+double __exp__D(x, c)
+double x, c;
+{
+       double  z,hi,lo, t;
+       int k;
+
+#if !defined(vax)&&!defined(tahoe)
+       if (x!=x) return(x);    /* x is NaN */
+#endif /* !defined(vax)&&!defined(tahoe) */
+       if ( x <= lnhuge ) {
+               if ( x >= lntiny ) {
+
+                   /* argument reduction : x --> x - k*ln2 */
+                       z = invln2*x;
+                       k = z + copysign(.5, x);
+
+                   /* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */
+
+                       hi=(x-k*ln2hi);                 /* Exact. */
+                       x= hi - (lo = k*ln2lo-c);
+                   /* return 2^k*[1+x+x*c/(2+c)]  */
+                       z=x*x;
+                       c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
+                       c = (x*c)/(2.0-c);
+
+                       return  scalb(1.+(hi-(lo - c)), k);
+               }
+               /* end of x > lntiny */
+
+               else 
+                    /* exp(-big#) underflows to zero */
+                    if(finite(x))  return(scalb(1.0,-5000));
+
+                    /* exp(-INF) is zero */
+                    else return(0.0);
+       }
+       /* end of x < lnhuge */
+
+       else 
+       /* exp(INF) is INF, exp(+big#) overflows to INF */
+           return( finite(x) ?  scalb(1.0,5000)  : x);
+}