projects / unix-history / blobdiff
? search:
summary | tags | clone url | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
BSD 4_3_Tahoe release
[unix-history] / usr / src / usr.lib / libm / exp.c
diff --git a/usr/src/usr.lib/libm/exp.c b/usr/src/usr.lib/libm/exp.c
index 701b2b5..242a2e2 100644 (file)
--- a/usr/src/usr.lib/libm/exp.c
+++ b/usr/src/usr.lib/libm/exp.c
@@ -1,45 +1,53 @@
-/* 
+/*
  * Copyright (c) 1985 Regents of the University of California.
  * Copyright (c) 1985 Regents of the University of California.
- * 
- * Use and reproduction of this software are granted  in  accordance  with
- * the terms and conditions specified in  the  Berkeley  Software  License
- * Agreement (in particular, this entails acknowledgement of the programs'
- * source, and inclusion of this notice) with the additional understanding
- * that  all  recipients  should regard themselves as participants  in  an
- * ongoing  research  project and hence should  feel  obligated  to report
- * their  experiences (good or bad) with these elementary function  codes,
- * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms are permitted
+ * provided that the above copyright notice and this paragraph are
+ * duplicated in all such forms and that any documentation,
+ * advertising materials, and other materials related to such
+ * distribution and use acknowledge that the software was developed
+ * by the University of California, Berkeley.  The name of the
+ * University may not be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
+ * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ * All recipients should regard themselves as participants in an ongoing
+ * research project and hence should feel obligated to report their
+ * experiences (good or bad) with these elementary function codes, using
+ * the sendbug(8) program, to the authors.
  */
 
 #ifndef lint
  */
 
 #ifndef lint
-static char sccsid[] = "@(#)exp.c      4.3 (Berkeley) 8/21/85";
-#endif not lint
+static char sccsid[] = "@(#)exp.c      5.3 (Berkeley) 6/30/88";
+#endif /* not lint */
 
 /* EXP(X)
  * RETURN THE EXPONENTIAL OF X
  * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
  * CODED IN C BY K.C. NG, 1/19/85; 
 
 /* EXP(X)
  * RETURN THE EXPONENTIAL OF X
  * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
  * CODED IN C BY K.C. NG, 1/19/85; 
- * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85.
+ * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
  *
  * Required system supported functions:
  *     scalb(x,n)      
  *     copysign(x,y)   
  *     finite(x)
  *
  *
  * Required system supported functions:
  *     scalb(x,n)      
  *     copysign(x,y)   
  *     finite(x)
  *
- * Kernel function:
- *     exp__E(x,c)
- *
  * Method:
  *     1. Argument Reduction: given the input x, find r and integer k such 
  *        that
  *                        x = k*ln2 + r,  |r| <= 0.5*ln2 .  
  *        r will be represented as r := z+c for better accuracy.
  *
  * Method:
  *     1. Argument Reduction: given the input x, find r and integer k such 
  *        that
  *                        x = k*ln2 + r,  |r| <= 0.5*ln2 .  
  *        r will be represented as r := z+c for better accuracy.
  *
- *     2. Compute expm1(r)=exp(r)-1 by 
+ *     2. Compute exp(r) by 
  *
  *
- *                     expm1(r=z+c) := z + exp__E(z,r)
+ *             exp(r) = 1 + r + r*R1/(2-R1),
+ *        where
+ *             R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))).
  *
  *
- *     3. exp(x) = 2^k * ( expm1(r) + 1 ).
+ *     3. exp(x) = 2^k * exp(r) .
  *
  * Special cases:
  *     exp(INF) is INF, exp(NaN) is NaN;
  *
  * Special cases:
  *     exp(INF) is INF, exp(NaN) is NaN;
@@ -49,7 +57,7 @@ static char sccsid[] = "@(#)exp.c     4.3 (Berkeley) 8/21/85";
  * Accuracy:
  *     exp(x) returns the exponential of x nearly rounded. In a test run
  *     with 1,156,000 random arguments on a VAX, the maximum observed
  * Accuracy:
  *     exp(x) returns the exponential of x nearly rounded. In a test run
  *     with 1,156,000 random arguments on a VAX, the maximum observed
- *     error was .768 ulps (units in the last place).
+ *     error was 0.869 ulps (units in the last place).
  *
  * Constants:
  * The hexadecimal values are the intended ones for the following constants.
  *
  * Constants:
  * The hexadecimal values are the intended ones for the following constants.
@@ -58,41 +66,67 @@ static char sccsid[] = "@(#)exp.c   4.3 (Berkeley) 8/21/85";
  * shown.
  */
 
  * shown.
  */
 
-#ifdef VAX     /* VAX D format */
-/* double static */
+#if defined(vax)||defined(tahoe)       /* VAX D format */
+#ifdef vax
+#define _0x(A,B)       0x/**/A/**/B
+#else  /* vax */
+#define _0x(A,B)       0x/**/B/**/A
+#endif /* vax */
+/* static double */
 /* ln2hi  =  6.9314718055829871446E-1    , Hex  2^  0   *  .B17217F7D00000 */
 /* ln2lo  =  1.6465949582897081279E-12   , Hex  2^-39   *  .E7BCD5E4F1D9CC */
 /* lnhuge =  9.4961163736712506989E1     , Hex  2^  7   *  .BDEC1DA73E9010 */
 /* lntiny = -9.5654310917272452386E1     , Hex  2^  7   * -.BF4F01D72E33AF */
 /* invln2 =  1.4426950408889634148E0     ; Hex  2^  1   *  .B8AA3B295C17F1 */
 /* ln2hi  =  6.9314718055829871446E-1    , Hex  2^  0   *  .B17217F7D00000 */
 /* ln2lo  =  1.6465949582897081279E-12   , Hex  2^-39   *  .E7BCD5E4F1D9CC */
 /* lnhuge =  9.4961163736712506989E1     , Hex  2^  7   *  .BDEC1DA73E9010 */
 /* lntiny = -9.5654310917272452386E1     , Hex  2^  7   * -.BF4F01D72E33AF */
 /* invln2 =  1.4426950408889634148E0     ; Hex  2^  1   *  .B8AA3B295C17F1 */
-static long     ln2hix[] = { 0x72174031, 0x0000f7d0};
-static long     ln2lox[] = { 0xbcd52ce7, 0xd9cce4f1};
-static long    lnhugex[] = { 0xec1d43bd, 0x9010a73e};
-static long    lntinyx[] = { 0x4f01c3bf, 0x33afd72e};
-static long    invln2x[] = { 0xaa3b40b8, 0x17f1295c};
+/* p1     =  1.6666666666666602251E-1    , Hex  2^-2    *  .AAAAAAAAAAA9F1 */
+/* p2     = -2.7777777777015591216E-3    , Hex  2^-8    * -.B60B60B5F5EC94 */
+/* p3     =  6.6137563214379341918E-5    , Hex  2^-13   *  .8AB355792EF15F */
+/* p4     = -1.6533902205465250480E-6    , Hex  2^-19   * -.DDEA0E2E935F84 */
+/* p5     =  4.1381367970572387085E-8    , Hex  2^-24   *  .B1BB4B95F52683 */
+static long     ln2hix[] = { _0x(7217,4031), _0x(0000,f7d0)};
+static long     ln2lox[] = { _0x(bcd5,2ce7), _0x(d9cc,e4f1)};
+static long    lnhugex[] = { _0x(ec1d,43bd), _0x(9010,a73e)};
+static long    lntinyx[] = { _0x(4f01,c3bf), _0x(33af,d72e)};
+static long    invln2x[] = { _0x(aa3b,40b8), _0x(17f1,295c)};
+static long        p1x[] = { _0x(aaaa,3f2a), _0x(a9f1,aaaa)};
+static long        p2x[] = { _0x(0b60,bc36), _0x(ec94,b5f5)};
+static long        p3x[] = { _0x(b355,398a), _0x(f15f,792e)};
+static long        p4x[] = { _0x(ea0e,b6dd), _0x(5f84,2e93)};
+static long        p5x[] = { _0x(bb4b,3431), _0x(2683,95f5)};
 #define    ln2hi    (*(double*)ln2hix)
 #define    ln2lo    (*(double*)ln2lox)
 #define   lnhuge    (*(double*)lnhugex)
 #define   lntiny    (*(double*)lntinyx)
 #define   invln2    (*(double*)invln2x)
 #define    ln2hi    (*(double*)ln2hix)
 #define    ln2lo    (*(double*)ln2lox)
 #define   lnhuge    (*(double*)lnhugex)
 #define   lntiny    (*(double*)lntinyx)
 #define   invln2    (*(double*)invln2x)
-#else  /* IEEE double */
-double static
+#define       p1    (*(double*)p1x)
+#define       p2    (*(double*)p2x)
+#define       p3    (*(double*)p3x)
+#define       p4    (*(double*)p4x)
+#define       p5    (*(double*)p5x)
+
+#else  /* defined(vax)||defined(tahoe) */
+static double
+p1     =  1.6666666666666601904E-1    , /*Hex  2^-3    *  1.555555555553E */
+p2     = -2.7777777777015593384E-3    , /*Hex  2^-9    * -1.6C16C16BEBD93 */
+p3     =  6.6137563214379343612E-5    , /*Hex  2^-14   *  1.1566AAF25DE2C */
+p4     = -1.6533902205465251539E-6    , /*Hex  2^-20   * -1.BBD41C5D26BF1 */
+p5     =  4.1381367970572384604E-8    , /*Hex  2^-25   *  1.6376972BEA4D0 */
 ln2hi  =  6.9314718036912381649E-1    , /*Hex  2^ -1   *  1.62E42FEE00000 */
 ln2lo  =  1.9082149292705877000E-10   , /*Hex  2^-33   *  1.A39EF35793C76 */
 lnhuge =  7.1602103751842355450E2     , /*Hex  2^  9   *  1.6602B15B7ECF2 */
 lntiny = -7.5137154372698068983E2     , /*Hex  2^  9   * -1.77AF8EBEAE354 */
 invln2 =  1.4426950408889633870E0     ; /*Hex  2^  0   *  1.71547652B82FE */
 ln2hi  =  6.9314718036912381649E-1    , /*Hex  2^ -1   *  1.62E42FEE00000 */
 ln2lo  =  1.9082149292705877000E-10   , /*Hex  2^-33   *  1.A39EF35793C76 */
 lnhuge =  7.1602103751842355450E2     , /*Hex  2^  9   *  1.6602B15B7ECF2 */
 lntiny = -7.5137154372698068983E2     , /*Hex  2^  9   * -1.77AF8EBEAE354 */
 invln2 =  1.4426950408889633870E0     ; /*Hex  2^  0   *  1.71547652B82FE */
-#endif
+#endif /* defined(vax)||defined(tahoe) */
 
 double exp(x)
 double x;
 {
 
 double exp(x)
 double x;
 {
-       double scalb(), copysign(), exp__E(), z,hi,lo,c;
+       double scalb(), copysign(), z,hi,lo,c;
        int k,finite();
 
        int k,finite();
 
-#ifndef VAX
+#if !defined(vax)&&!defined(tahoe)
        if(x!=x) return(x);     /* x is NaN */
        if(x!=x) return(x);     /* x is NaN */
-#endif
+#endif /* !defined(vax)&&!defined(tahoe) */
        if( x <= lnhuge ) {
                if( x >= lntiny ) {
 
        if( x <= lnhuge ) {
                if( x >= lntiny ) {
 
@@ -100,14 +134,16 @@ double x;
 
                        k=invln2*x+copysign(0.5,x);     /* k=NINT(x/ln2) */
 
 
                        k=invln2*x+copysign(0.5,x);     /* k=NINT(x/ln2) */
 
-                       /* express x-k*ln2 as z+c */
+                   /* express x-k*ln2 as hi-lo and let x=hi-lo rounded */
+
                        hi=x-k*ln2hi;
                        hi=x-k*ln2hi;
-                       z=hi-(lo=k*ln2lo);
-                       c=(hi-z)-lo;
+                       x=hi-(lo=k*ln2lo);
+
+                   /* return 2^k*[1+x+x*c/(2+c)]  */
+                       z=x*x;
+                       c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
+                       return  scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k);
 
 
-                   /* return 2^k*[expm1(x) + 1]  */
-                       z += exp__E(z,c);
-                       return (scalb(z+1.0,k));  
                }
                /* end of x > lntiny */
 
                }
                /* end of x > lntiny */
 
Continuous source code history from original PDP-7 UNIX to FreeBSD 2, the first fully unencumbered FreeBSD release.