OpenE2K
/
glibc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Optimized version of powf() 

Most significant changes are code simplification and use of doubles for
intermediate values.  Also, some rearrangement to move early
non-dependent code later, out of the faster paths.

	* sysdeps/ieee754/flt-32/e_powf.c: Optimized implementation utilizing
	rearranged code and doubles float types.
This commit is contained in:
Paul Clarke 2017-06-23 09:10:32 -03:00 committed by Tulio Magno Quites Machado Filho
parent 76c4c838e7
commit c340290de8
2 changed files with 75 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
ChangeLog
View File

@ -1,3 +1,8 @@
2017-06-23 Paul A. Clarke <pc@us.ibm.com>
* sysdeps/ieee754/flt-32/e_powf.c: Optimized implementation utilizing
rearranged code and doubles float types.
2017-06-23 Rogerio A. Cardoso <rcardoso@linux.vnet.ibm.com>
* sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-ppc64.c: Change

171
sysdeps/ieee754/flt-32/e_powf.c
View File

@ -1,6 +1,22 @@
/* e_powf.c -- float version of e_pow.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/* Copyright (C) 2017 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<http://www.gnu.org/licenses/>. */
/*
* ====================================================
@ -20,8 +36,6 @@ 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 */
dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
zero = 0.0,
one = 1.0,
two = 2.0,
@ -38,28 +52,24 @@ P2 = -2.7777778450e-03, /* 0xbb360b61 */
P3 = 6.6137559770e-05, /* 0x388ab355 */
P4 = -1.6533901999e-06, /* 0xb5ddea0e */
P5 = 4.1381369442e-08, /* 0x3331bb4c */
lg2 = 6.9314718246e-01, /* 0x3f317218 */
lg2_h = 6.93145752e-01, /* 0x3f317200 */
lg2_l = 1.42860654e-06, /* 0x35bfbe8c */
ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
cp_h = 0xf.64p-4, /* cp high 12 bits. */
cp_l = -0x7.b11e3p-16, /* 2/(3ln2) - 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*/
ovt = 4.2995665694e-08; /* -(128-log2(ovfl+.5ulp)) */
static const double
dp[] = { 0.0, 0x1.2b803473f7ad1p-1, }, /* log2(1.5) */
lg2 = M_LN2,
cp = 2.0/3.0/M_LN2,
invln2 = 1.0/M_LN2;
float
__ieee754_powf(float x, float y)
{
float z,ax,z_h,z_l,p_h,p_l;
float y1,t1,t2,r,s,t,u,v,w;
float z, ax, s;
double d1, d2;
int32_t i,j,k,yisint,n;
int32_t hx,hy,ix,iy,is;
int32_t hx,hy,ix,iy;
GET_FLOAT_WORD(hx,x);
GET_FLOAT_WORD(hy,y);
ix = hx&0x7fffffff; iy = hy&0x7fffffff;
iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
if(iy==0 && !issignaling (x)) return one;
@ -68,26 +78,14 @@ __ieee754_powf(float x, float y)
if(x == 1.0 && !issignaling (y)) return one;
if(x == -1.0 && isinf(y)) return one;
GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff;
/* +-NaN return x+y */
if(__builtin_expect(ix > 0x7f800000 ||
iy > 0x7f800000, 0))
return x+y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
yisint = 0;
if(hx<0) {
if(iy>=0x4b800000) yisint = 2; /* even integer y */
else if(iy>=0x3f800000) {
k = (iy>>23)-0x7f; /* exponent */
j = iy>>(23-k);
if((j<<(23-k))==iy) yisint = 2-(j&1);
}
}
/* special value of y */
if (__builtin_expect(iy==0x7f800000, 0)) { /* y is +-inf */
if (ix==0x3f800000)
@ -106,6 +104,21 @@ __ieee754_powf(float x, float y)
return __ieee754_sqrtf(x);
}
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
yisint = 0;
if(hx<0) {
if(iy>=0x4b800000) yisint = 2; /* even integer y */
else if(iy>=0x3f800000) {
k = (iy>>23)-0x7f; /* exponent */
j = iy>>(23-k);
if((j<<(23-k))==iy) yisint = 2-(j&1);
}
}
ax = fabsf(x);
/* special value of x */
if(__builtin_expect(ix==0x7f800000||ix==0||ix==0x3f800000, 0)){
@ -131,16 +144,10 @@ __ieee754_powf(float x, float y)
if(ix>0x3f800007) return (hy>0)? huge*huge: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-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;
t1 = u+v;
GET_FLOAT_WORD(is,t1);
SET_FLOAT_WORD(t1,is&0xfffff000);
t2 = v-(t1-u);
d2 = ax-1; /* d2 has 20 trailing zeros. */
d2 = d2 * invln2 -
(d2 * d2) * (0.5 - d2 * (0.333333333333 - d2 * 0.25)) * invln2;
} else {
float s2,s_h,s_l,t_h,t_l;
/* Avoid internal underflow for tiny y. The exact value
of y does not matter if |y| <= 2**-32. */
if (iy < 0x2f800000)
@ -158,69 +165,37 @@ __ieee754_powf(float x, float y)
else {k=0;n+=1;ix -= 0x00800000;}
SET_FLOAT_WORD(ax,ix);
/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]);
s = u*v;
s_h = s;
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)+0x00400000+(k<<21))
& 0xfffff000));
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
/* compute log(ax) */
s2 = s*s;
r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
r += s_l*(s_h+s);
s2 = s_h*s_h;
t_h = (float)3.0+s2+r;
GET_FLOAT_WORD(is,t_h);
SET_FLOAT_WORD(t_h,is&0xfffff000);
t_l = r-((t_h-(float)3.0)-s2);
/* u+v = s*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*s;
/* 2/(3log2)*(s+...) */
p_h = u+v;
GET_FLOAT_WORD(is,p_h);
SET_FLOAT_WORD(p_h,is&0xfffff000);
p_l = v-(p_h-u);
z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l*p_h+p_l*cp+dp_l[k];
/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
t = (float)n;
t1 = (((z_h+z_l)+dp_h[k])+t);
GET_FLOAT_WORD(is,t1);
SET_FLOAT_WORD(t1,is&0xfffff000);
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
/* compute d1 = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
d1 = (ax-(double)bp[k])/(ax+(double)bp[k]);
/* compute d2 = log(ax) */
d2 = d1 * d1;
d2 = 3.0 + d2 + d2*d2*(L1+d2*(L2+d2*(L3+d2*(L4+d2*(L5+d2*L6)))));
/* 2/(3log2)*(d2+...) */
d2 = d1*d2*cp;
/* log2(ax) = (d2+..)*2/(3*log2) */
d2 = d2+dp[k]+(double)n;
}
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 y1+y2 and compute (y1+y2)*(t1+t2) */
GET_FLOAT_WORD(is,y);
SET_FLOAT_WORD(y1,is&0xfffff000);
p_l = (y-y1)*t1+y*t2;
p_h = y1*t1;
z = p_l+p_h;
/* compute y * d2 */
d1 = y * d2;
z = d1;
GET_FLOAT_WORD(j,z);
if (__builtin_expect(j>0x43000000, 0)) /* if z > 128 */
return s*huge*huge; /* overflow */
else if (__builtin_expect(j==0x43000000, 0)) { /* if z == 128 */
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
if(ovt>(z-d1)) return s*huge*huge; /* overflow */
}
else if (__builtin_expect((j&0x7fffffff)>0x43160000, 0))/* z <= -150 */
return s*tiny*tiny; /* underflow */
else if (__builtin_expect((u_int32_t) j==0xc3160000, 0)){/* z == -150*/
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
if(0.0<=(z-d1)) return s*tiny*tiny; /* underflow */
}
/*
* compute 2**(p_h+p_l)
* compute 2**d1
*/
i = j&0x7fffffff;
k = (i>>23)-0x7f;
@ -228,22 +203,16 @@ __ieee754_powf(float x, float y)
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 */
SET_FLOAT_WORD(t,n&~(0x007fffff>>k));
SET_FLOAT_WORD(z,n&~(0x007fffff>>k));
n = ((n&0x007fffff)|0x00800000)>>(23-k);
if(j<0) n = -n;
p_h -= t;
d1 -= z;
}
t = p_l+p_h;
GET_FLOAT_WORD(is,t);
SET_FLOAT_WORD(t,is&0xfffff000);
u = t*lg2_h;
v = (p_l-(t-p_h))*lg2+t*lg2_l;
z = u+v;
w = v-(z-u);
t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two)-(w+z*w);
z = one-(r-z);
d1 = d1 * lg2;
d2 = d1*d1;
d2 = d1 - d2*(P1+d2*(P2+d2*(P3+d2*(P4+d2*P5))));
d2 = (d1*d2)/(d2-two);
z = one - (d2-d1);
GET_FLOAT_WORD(j,z);
j += (n<<23);
if((j>>23)<=0) /* subnormal output */