From 6d3144f9382a4aede51913d7465d300a320412e3 Mon Sep 17 00:00:00 2001 From: Raymond Toy <toy.raymond@gmail.com> Date: Fri, 1 Aug 2014 19:54:27 -0700 Subject: [PATCH] Update to use unions. * e_pow.c: * Use unions to get to the parts of a double * fdlibm.h: * Declare sqrt function. --- src/lisp/e_pow.c | 71 +++++++++++++++++++++++++++++++++++------------ src/lisp/fdlibm.h | 1 + 2 files changed, 54 insertions(+), 18 deletions(-) diff --git a/src/lisp/e_pow.c b/src/lisp/e_pow.c index 5683bf5fd..9dafc2649 100644 --- a/src/lisp/e_pow.c +++ b/src/lisp/e_pow.c @@ -109,10 +109,15 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ int i0,i1,i,j,k,yisint,n; int hx,hy,ix,iy; unsigned lx,ly; + union { int i[2]; double d; } ux; + union { int i[2]; double d; } uy; + union { int i[2]; double d; } utmp; i0 = ((*(int*)&one)>>29)^1; i1=1-i0; - hx = __HI(x); lx = __LO(x); - hy = __HI(y); ly = __LO(y); + ux.d = x; + hx = ux.i[HIWORD]; lx = ux.i[LOWORD]; + uy.d = y; + hy = uy.i[HIWORD]; ly = uy.i[LOWORD]; ix = hx&0x7fffffff; iy = hy&0x7fffffff; /* y==zero: x**0 = 1 */ @@ -203,14 +208,20 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ v = t*ivln2_l-w*ivln2; t1 = u+v; - __LO(t1) = 0; + utmp.d = t1; + utmp.i[LOWORD] = 0; + t1 = utmp.d; t2 = v-(t1-u); } else { double ss,s2,s_h,s_l,t_h,t_l; n = 0; /* take care subnormal number */ - if(ix<0x00100000) - {ax *= two53; n -= 53; ix = __HI(ax); } + if(ix<0x00100000) { + ax *= two53; + n -= 53; + utmp.d = ax; + ix = utmp.i[HIWORD]; + } n += ((ix)>>20)-0x3ff; j = ix&0x000fffff; /* determine interval */ @@ -218,17 +229,23 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ else {k=0;n+=1;ix -= 0x00100000;} - __HI(ax) = ix; + utmp.d = ax; + utmp.i[HIWORD] = ix; + ax = utmp.d; /* compute ss = 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]); ss = u*v; s_h = ss; - __LO(s_h) = 0; + utmp.d = s_h; + utmp.i[LOWORD] = 0; + s_h = utmp.d; /* t_h=ax+bp[k] High */ t_h = zero; - __HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18); + utmp.d = t_h; + utmp.i[HIWORD]=((ix>>1)|0x20000000)+0x00080000+(k<<18); + t_h = utmp.d; t_l = ax - (t_h-bp[k]); s_l = v*((u-s_h*t_h)-s_h*t_l); /* compute log(ax) */ @@ -237,32 +254,41 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ r += s_l*(s_h+ss); s2 = s_h*s_h; t_h = 3.0+s2+r; - __LO(t_h) = 0; + utmp.d = t_h; + utmp.i[LOWORD] = 0; + t_h = utmp.d; t_l = r-((t_h-3.0)-s2); /* u+v = ss*(1+...) */ u = s_h*t_h; v = s_l*t_h+t_l*ss; /* 2/(3log2)*(ss+...) */ p_h = u+v; - __LO(p_h) = 0; + utmp.d = p_h; + utmp.i[LOWORD] = 0; + p_h = utmp.d; 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) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ t = (double)n; t1 = (((z_h+z_l)+dp_h[k])+t); - __LO(t1) = 0; + utmp.d = t1; + utmp.i[LOWORD] = 0; + t1 = utmp.d; t2 = z_l-(((t1-t)-dp_h[k])-z_h); } /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ y1 = y; - __LO(y1) = 0; + utmp.d = y1; + utmp.i[LOWORD] = 0; + y1 = utmp.d; p_l = (y-y1)*t1+y*t2; p_h = y1*t1; z = p_l+p_h; - j = __HI(z); - i = __LO(z); + utmp.d = z; + j = utmp.i[HIWORD]; + i = utmp.i[LOWORD]; if (j>=0x40900000) { /* z >= 1024 */ if(((j-0x40900000)|i)!=0) /* if z > 1024 */ return s*huge*huge; /* overflow */ @@ -286,13 +312,17 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ n = j+(0x00100000>>(k+1)); k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */ t = zero; - __HI(t) = (n&~(0x000fffff>>k)); + utmp.d = t; + utmp.i[HIWORD] = (n&~(0x000fffff>>k)); + t = utmp.d; n = ((n&0x000fffff)|0x00100000)>>(20-k); if(j<0) n = -n; p_h -= t; } t = p_l+p_h; - __LO(t) = 0; + utmp.d = t; + utmp.i[LOWORD] = 0; + t = utmp.d; u = t*lg2_h; v = (p_l-(t-p_h))*lg2+t*lg2_l; z = u+v; @@ -301,9 +331,14 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); r = (z*t1)/(t1-two)-(w+z*w); z = one-(r-z); - j = __HI(z); + utmp.d = z; + j = utmp.i[HIWORD]; j += (n<<20); if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */ - else __HI(z) += (n<<20); + else { + utmp.d = z; + utmp.i[HIWORD] += (n<<20); + z = utmp.d; + } return s*z; } diff --git a/src/lisp/fdlibm.h b/src/lisp/fdlibm.h index 467daa537..180137baf 100644 --- a/src/lisp/fdlibm.h +++ b/src/lisp/fdlibm.h @@ -32,6 +32,7 @@ enum { HIWORD = 0, LOWORD = 1 }; */ extern double fabs(double); extern double floor(double); +extern double sqrt(double); /* ieee style elementary functions */ extern int __ieee754_rem_pio2(double,double*); -- GitLab