BSD 4_4 release

[unix-history] / usr / src / lib / libm / common_source / exp.c

diff --git a/usr/src/lib/libm/common_source/exp.c b/usr/src/lib/libm/common_source/exp.c
index cd97eff..9b4f045 100644 (file)
--- a/usr/src/lib/libm/common_source/exp.c
+++ b/usr/src/lib/libm/common_source/exp.c
@@ -1,12 +1,38 @@
 /*
 /*

- * 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.

  *
  *

- * %sccs.include.redist.c%
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. All advertising materials mentioning features or use of this software
+ *    must display the following acknowledgement:
+ *     This product includes software developed by the University of
+ *     California, Berkeley and its contributors.
+ * 4. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.

  */
 
 #ifndef lint
  */
 
 #ifndef lint

-static char sccsid[] = "@(#)exp.c      5.6 (Berkeley) %G%";
+static char sccsid[] = "@(#)exp.c      8.1 (Berkeley) 6/4/93";

 #endif /* not lint */
 
 /* EXP(X)
 #endif /* not lint */
 
 /* EXP(X)


@@ -130,3 +156,48 @@ double x;
        /* exp(INF) is INF, exp(+big#) overflows to INF */
            return( finite(x) ?  scalb(1.0,5000)  : 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);
+}