Source-Changes-HG archive

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index][Old Index]

[src/trunk]: src/lib/libm/src sync with FreeBSD



details:   https://anonhg.NetBSD.org/src/rev/03f23d43ed6b
branches:  trunk
changeset: 817616:03f23d43ed6b
user:      christos <christos%NetBSD.org@localhost>
date:      Sat Aug 27 10:01:08 2016 +0000

description:
sync with FreeBSD

diffstat:

 lib/libm/src/e_pow.c  |  24 +++++++++++------
 lib/libm/src/e_powf.c |  69 +++++++++++++++++++++++++++-----------------------
 2 files changed, 52 insertions(+), 41 deletions(-)

diffs (269 lines):

diff -r 303a08c32450 -r 03f23d43ed6b lib/libm/src/e_pow.c
--- a/lib/libm/src/e_pow.c      Sat Aug 27 10:00:38 2016 +0000
+++ b/lib/libm/src/e_pow.c      Sat Aug 27 10:01:08 2016 +0000
@@ -1,9 +1,8 @@
-/* @(#)e_pow.c 5.1 93/09/24 */
+/* @(#)e_pow.c 1.5 04/04/22 SMI */
 /*
  * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
  * software is freely granted, provided that this notice
  * is preserved.
@@ -12,7 +11,7 @@
 
 #include <sys/cdefs.h>
 #if defined(LIBM_SCCS) && !defined(lint)
-__RCSID("$NetBSD: e_pow.c,v 1.16 2010/04/23 19:17:07 drochner Exp $");
+__RCSID("$NetBSD: e_pow.c,v 1.17 2016/08/27 10:01:08 christos Exp $");
 #endif
 
 /* __ieee754_pow(x,y) return x**y
@@ -29,13 +28,13 @@
  * Special cases:
  *     1.  (anything) ** 0  is 1
  *     2.  (anything) ** 1  is itself
- *     3.  (anything) ** NAN is NAN
+ *     3.  (anything) ** NAN is NAN except 1 ** NAN = 1
  *     4.  NAN ** (anything except 0) is NAN
  *     5.  +-(|x| > 1) **  +INF is +INF
  *     6.  +-(|x| > 1) **  -INF is +0
  *     7.  +-(|x| < 1) **  +INF is +0
  *     8.  +-(|x| < 1) **  -INF is +INF
- *     9.  +-1         ** +-INF is NAN
+ *     9.  +-1         ** +-INF is 1
  *     10. +0 ** (+anything except 0, NAN)               is +0
  *     11. -0 ** (+anything except 0, NAN, odd integer)  is +0
  *     12. +0 ** (-anything except 0, NAN)               is +INF
@@ -113,10 +112,13 @@
     /* y==zero: x**0 = 1 */
        if((iy|ly)==0) return one;
 
-    /* +-NaN return x+y */
+    /* x==1: 1**y = 1, even if y is NaN */
+       if (hx==0x3ff00000 && lx == 0) return one;
+
+    /* y!=zero: result is NaN if either arg is NaN */
        if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
           iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
-               return x+y;
+               return (x+0.0)+(y+0.0);
 
     /* determine if y is an odd int when x < 0
      * yisint = 0      ... y is not an integer
@@ -142,7 +144,7 @@
        if(ly==0) {
            if (iy==0x7ff00000) {       /* y is +-inf */
                if(((ix-0x3ff00000)|lx)==0)
-                   return  y - y;      /* inf**+-1 is NaN */
+                   return  one;        /* (-1)**+-inf is 1 */
                else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
                    return (hy>=0)? y: zero;
                else                    /* (|x|<1)**-,+inf = inf,0 */
@@ -174,7 +176,11 @@
            }
        }
 
+    /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be
        n = (hx>>31)+1;
+       but ANSI C says a right shift of a signed negative quantity is
+       implementation defined.  */
+       n = ((u_int32_t)hx>>31)-1;
 
     /* (x<0)**(non-int) is NaN */
        if((n|yisint)==0) return (x-x)/(x-x);
diff -r 303a08c32450 -r 03f23d43ed6b lib/libm/src/e_powf.c
--- a/lib/libm/src/e_powf.c     Sat Aug 27 10:00:38 2016 +0000
+++ b/lib/libm/src/e_powf.c     Sat Aug 27 10:01:08 2016 +0000
@@ -15,15 +15,13 @@
 
 #include <sys/cdefs.h>
 #if defined(LIBM_SCCS) && !defined(lint)
-__RCSID("$NetBSD: e_powf.c,v 1.15 2010/04/23 19:17:07 drochner Exp $");
+__RCSID("$NetBSD: e_powf.c,v 1.16 2016/08/27 10:01:08 christos Exp $");
 #endif
 
 #include "namespace.h"
 #include "math.h"
 #include "math_private.h"
 
-static const float huge = 1.0e+30, tiny = 1.0e-30;
-
 static const float
 bp[] = {1.0, 1.5,},
 dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
@@ -32,6 +30,8 @@
 one    =  1.0,
 two    =  2.0,
 two24  =  16777216.0,  /* 0x4b800000 */
+huge   =  1.0e30,
+tiny    =  1.0e-30,
        /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
 L1  =  6.0000002384e-01, /* 0x3f19999a */
 L2  =  4.2857143283e-01, /* 0x3edb6db7 */
@@ -49,8 +49,8 @@
 lg2_l  =  1.42860654e-06, /* 0x35bfbe8c */
 ovt =  4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
 cp    =  9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
-cp_h  =  9.6179199219e-01, /* 0x3f763800 =head of cp */
-cp_l  =  4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */
+cp_h  =  9.6191406250e-01, /* 0x3f764000 =12b cp */
+cp_l  = -1.1736857402e-04, /* 0xb8f623c6 =tail of cp_h */
 ivln2    =  1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
 ivln2_h  =  1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
 ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
@@ -59,7 +59,7 @@
 __ieee754_powf(float x, float y)
 {
        float z,ax,z_h,z_l,p_h,p_l;
-       float yy1,t1,t2,r,s,t,u,v,w;
+       float yy1,t1,t2,r,s,sn,t,u,v,w;
        int32_t i,j,k,yisint,n;
        int32_t hx,hy,ix,iy,is;
 
@@ -70,10 +70,13 @@
     /* y==zero: x**0 = 1 */
        if(iy==0) return one;
 
-    /* +-NaN return x+y */
+    /* x==1: 1**y = 1, even if y is NaN */
+       if (hx==0x3f800000) return one;
+
+    /* y!=zero: result is NaN if either arg is NaN */
        if(ix > 0x7f800000 ||
           iy > 0x7f800000)
-               return x+y;
+               return (x+0.0F)+(y+0.0F);
 
     /* determine if y is an odd int when x < 0
      * yisint = 0      ... y is not an integer
@@ -84,8 +87,8 @@
        if(hx<0) {
            if(iy>=0x4b800000) yisint = 2; /* even integer y */
            else if(iy>=0x3f800000) {
-               k = (iy>>23)-0x7f;         /* exponent */
-               j = iy>>(23-k);
+               k = ((uint32_t)iy>>23)-0x7f;       /* exponent */
+               j = (uint32_t)iy>>(23-k);
                if((j<<(23-k))==iy) yisint = 2-(j&1);
            }
        }
@@ -93,7 +96,7 @@
     /* special value of y */
        if (iy==0x7f800000) {   /* y is +-inf */
            if (ix==0x3f800000)
-               return  y - y;  /* inf**+-1 is NaN */
+               return  one;    /* (-1)**+-inf is NaN */
            else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
                return (hy>=0)? y: zero;
            else                        /* (|x|<1)**-,+inf = inf,0 */
@@ -122,17 +125,22 @@
            return z;
        }
 
+       n = ((u_int32_t)hx>>31)-1;
+
     /* (x<0)**(non-int) is NaN */
-       if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
+       if((n|yisint)==0) return (x-x)/(x-x);
+
+       sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+       if((n|(yisint-1))==0) sn = -one;/* (-ve)**(odd int) */
 
     /* |y| is huge */
        if(iy>0x4d000000) { /* if |y| > 2**27 */
        /* over/underflow if x is not close to one */
-           if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny;
-           if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny;
+           if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+           if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
        /* now |1-x| is tiny <= 2**-20, suffice to compute
           log(x) by x-x^2/2+x^3/3-x^4/4 */
-           t = ax-one;         /* t has 20 trailing zeros */
+           t = ax-1;           /* t has 20 trailing zeros */
            w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
            u = ivln2_h*t;      /* ivln2_h has 16 sig. bits */
            v = t*ivln2_l-w*ivln2;
@@ -146,7 +154,7 @@
        /* take care subnormal number */
            if(ix<0x00800000)
                {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); }
-           n  += ((ix)>>23)-0x7f;
+           n  += (((uint32_t)ix)>>23)-0x7f;
            j  = ix&0x007fffff;
        /* determine interval */
            ix = j|0x3f800000;          /* normalize ix */
@@ -163,7 +171,8 @@
            GET_FLOAT_WORD(is,s_h);
            SET_FLOAT_WORD(s_h,is&0xfffff000);
        /* t_h=ax+bp[k] High */
-           SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21));
+           is = (((uint32_t)ix>>1)&0xfffff000)|0x20000000;
+           SET_FLOAT_WORD(t_h,is+0x00400000+(k<<21));
            t_l = ax - (t_h-bp[k]);
            s_l = v*((u-s_h*t_h)-s_h*t_l);
        /* compute log(ax) */
@@ -193,10 +202,6 @@
            t2 = z_l-(((t1-t)-dp_h[k])-z_h);
        }
 
-       s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
-       if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
-           s = -one;   /* (-ve)**(odd int) */
-
     /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */
        GET_FLOAT_WORD(is,y);
        SET_FLOAT_WORD(yy1,is&0xfffff000);
@@ -205,32 +210,32 @@
        z = p_l+p_h;
        GET_FLOAT_WORD(j,z);
        if (j>0x43000000)                               /* if z > 128 */
-           return s*huge*huge;                         /* overflow */
+           return sn*huge*huge;                        /* overflow */
        else if (j==0x43000000) {                       /* if z == 128 */
-           if(p_l+ovt>z-p_h) return s*huge*huge;       /* overflow */
-       }
-       else if ((uint32_t)j==0xc3160000){              /* z == -150 */
-           if(p_l<=z-p_h) return s*tiny*tiny;          /* underflow */
+           if(p_l+ovt>z-p_h) return sn*huge*huge;      /* overflow */
        }
        else if ((j&0x7fffffff)>0x43160000)             /* z <= -150 */
-           return s*tiny*tiny;                         /* underflow */
+           return sn*tiny*tiny;                        /* underflow */
+       else if (j==(int)0xc3160000){                   /* z == -150 */
+           if(p_l<=z-p_h) return sn*tiny*tiny;         /* underflow */
+       }
     /*
      * compute 2**(p_h+p_l)
      */
        i = j&0x7fffffff;
-       k = (i>>23)-0x7f;
+       k = ((uint32_t)i>>23)-0x7f;
        n = 0;
        if(i>0x3f000000) {              /* if |z| > 0.5, set n = [z+0.5] */
            n = j+(0x00800000>>(k+1));
-           k = ((n&0x7fffffff)>>23)-0x7f;      /* new k for n */
+           k = (((uint32_t)n&0x7fffffff)>>23)-0x7f;    /* new k for n */
            SET_FLOAT_WORD(t,n&~(0x007fffff>>k));
-           n = ((n&0x007fffff)|0x00800000)>>(23-k);
+           n = (((uint32_t)n&0x007fffff)|0x00800000)>>(23-k);
            if(j<0) n = -n;
            p_h -= t;
        }
        t = p_l+p_h;
        GET_FLOAT_WORD(is,t);
-       SET_FLOAT_WORD(t,is&0xfffff000);
+       SET_FLOAT_WORD(t,is&0xffff8000);
        u = t*lg2_h;
        v = (p_l-(t-p_h))*lg2+t*lg2_l;
        z = u+v;
@@ -243,5 +248,5 @@
        j += (n<<23);
        if((j>>23)<=0) z = scalbnf(z,n);        /* subnormal output */
        else SET_FLOAT_WORD(z,j);
-       return s*z;
+       return sn*z;
 }



Home | Main Index | Thread Index | Old Index