diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/dla.h')
-rw-r--r-- | sysdeps/ieee754/dbl-64/dla.h | 119 |
1 files changed, 67 insertions, 52 deletions
diff --git a/sysdeps/ieee754/dbl-64/dla.h b/sysdeps/ieee754/dbl-64/dla.h index bf73fa902e..9f095f9bf5 100644 --- a/sysdeps/ieee754/dbl-64/dla.h +++ b/sysdeps/ieee754/dbl-64/dla.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -35,6 +35,18 @@ /* IEEE double. */ /***********************************************************************/ +/* We can use fma instructions if available. */ +#if defined __x86_64__ || (defined __i386__ && defined __SSE2_MATH__) +# ifdef __FMA4__ +# define DLA_FMA(x,y,z) \ + ({ double __zz; \ + asm ("vfmaddsd %3, %2, %1, %0" \ + : "=x" (__zz) : "x" (x), "xm" (y), "x" (-z)); \ + __zz; }) +# endif +#endif + + /* CN = 1+2**27 = '41a0000002000000' IEEE double format */ #define CN 134217729.0 @@ -44,7 +56,7 @@ /* z+zz = x+y exactly. */ #define EADD(x,y,z,zz) \ - z=(x)+(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x)); + z=(x)+(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x)); /* Exact subtraction of two single-length floating point numbers, Dekker. */ @@ -52,7 +64,7 @@ /* z+zz = x-y exactly. */ #define ESUB(x,y,z,zz) \ - z=(x)-(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); + z=(x)-(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); /* Exact multiplication of two single-length floating point numbers, */ @@ -60,10 +72,15 @@ /* satisfies z+zz = x*y exactly. p,hx,tx,hy,ty are temporary */ /* storage variables of type double. */ -#define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \ - p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \ - p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \ - z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty; +#ifdef DLA_FMA +# define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \ + z=x*y; zz=DLA_FMA(x,y,z); +#else +# define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \ + p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \ + p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \ + z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty; +#endif /* Exact multiplication of two single-length floating point numbers, Dekker. */ @@ -71,10 +88,15 @@ /* that satisfies z+zz = x*y exactly. p,hx,tx,hy,ty,q are temporary */ /* storage variables of type double. */ -#define MUL12(x,y,z,zz,p,hx,tx,hy,ty,q) \ - p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \ - p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \ - p=hx*hy; q=hx*ty+tx*hy; z=p+q; zz=((p-z)+q)+tx*ty; +#ifdef DLA_FMA +# define MUL12(x,y,z,zz,p,hx,tx,hy,ty,q) \ + EMULV(x,y,z,zz,p,hx,tx,hy,ty) +#else +# define MUL12(x,y,z,zz,p,hx,tx,hy,ty,q) \ + p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \ + p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \ + p=hx*hy; q=hx*ty+tx*hy; z=p+q; zz=((p-z)+q)+tx*ty; +#endif /* Double-length addition, Dekker. The macro produces a double-length */ @@ -84,10 +106,10 @@ /* storage variables of type double. */ #define ADD2(x,xx,y,yy,z,zz,r,s) \ - r=(x)+(y); s=(ABS(x)>ABS(y)) ? \ - (((((x)-r)+(y))+(yy))+(xx)) : \ - (((((y)-r)+(x))+(xx))+(yy)); \ - z=r+s; zz=(r-z)+s; + r=(x)+(y); s=(ABS(x)>ABS(y)) ? \ + (((((x)-r)+(y))+(yy))+(xx)) : \ + (((((y)-r)+(x))+(xx))+(yy)); \ + z=r+s; zz=(r-z)+s; /* Double-length subtraction, Dekker. The macro produces a double-length */ @@ -97,10 +119,10 @@ /* storage variables of type double. */ #define SUB2(x,xx,y,yy,z,zz,r,s) \ - r=(x)-(y); s=(ABS(x)>ABS(y)) ? \ - (((((x)-r)-(y))-(yy))+(xx)) : \ - ((((x)-((y)+r))+(xx))-(yy)); \ - z=r+s; zz=(r-z)+s; + r=(x)-(y); s=(ABS(x)>ABS(y)) ? \ + (((((x)-r)-(y))-(yy))+(xx)) : \ + ((((x)-((y)+r))+(xx))-(yy)); \ + z=r+s; zz=(r-z)+s; /* Double-length multiplication, Dekker. The macro produces a double-length */ @@ -110,8 +132,8 @@ /* temporary storage variables of type double. */ #define MUL2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc) \ - MUL12(x,y,c,cc,p,hx,tx,hy,ty,q) \ - cc=((x)*(yy)+(xx)*(y))+cc; z=c+cc; zz=(c-z)+cc; + MUL12(x,y,c,cc,p,hx,tx,hy,ty,q) \ + cc=((x)*(yy)+(xx)*(y))+cc; z=c+cc; zz=(c-z)+cc; /* Double-length division, Dekker. The macro produces a double-length */ @@ -121,8 +143,8 @@ /* are temporary storage variables of type double. */ #define DIV2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc,u,uu) \ - c=(x)/(y); MUL12(c,y,u,uu,p,hx,tx,hy,ty,q) \ - cc=(((((x)-u)-uu)+(xx))-c*(yy))/(y); z=c+cc; zz=(c-z)+cc; + c=(x)/(y); MUL12(c,y,u,uu,p,hx,tx,hy,ty,q) \ + cc=(((((x)-u)-uu)+(xx))-c*(yy))/(y); z=c+cc; zz=(c-z)+cc; /* Double-length addition, slower but more accurate than ADD2. */ @@ -133,17 +155,17 @@ /* are temporary storage variables of type double. */ #define ADD2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \ - r=(x)+(y); \ - if (ABS(x)>ABS(y)) { rr=((x)-r)+(y); s=(rr+(yy))+(xx); } \ - else { rr=((y)-r)+(x); s=(rr+(xx))+(yy); } \ - if (rr!=0.0) { \ - z=r+s; zz=(r-z)+s; } \ - else { \ - ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \ - u=r+s; \ - uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \ - w=uu+ss; z=u+w; \ - zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; } + r=(x)+(y); \ + if (ABS(x)>ABS(y)) { rr=((x)-r)+(y); s=(rr+(yy))+(xx); } \ + else { rr=((y)-r)+(x); s=(rr+(xx))+(yy); } \ + if (rr!=0.0) { \ + z=r+s; zz=(r-z)+s; } \ + else { \ + ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \ + u=r+s; \ + uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \ + w=uu+ss; z=u+w; \ + zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; } /* Double-length subtraction, slower but more accurate than SUB2. */ @@ -154,21 +176,14 @@ /* are temporary storage variables of type double. */ #define SUB2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \ - r=(x)-(y); \ - if (ABS(x)>ABS(y)) { rr=((x)-r)-(y); s=(rr-(yy))+(xx); } \ - else { rr=(x)-((y)+r); s=(rr+(xx))-(yy); } \ - if (rr!=0.0) { \ - z=r+s; zz=(r-z)+s; } \ - else { \ - ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \ - u=r+s; \ - uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \ - w=uu+ss; z=u+w; \ - zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; } - - - - - - - + r=(x)-(y); \ + if (ABS(x)>ABS(y)) { rr=((x)-r)-(y); s=(rr-(yy))+(xx); } \ + else { rr=(x)-((y)+r); s=(rr+(xx))-(yy); } \ + if (rr!=0.0) { \ + z=r+s; zz=(r-z)+s; } \ + else { \ + ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \ + u=r+s; \ + uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \ + w=uu+ss; z=u+w; \ + zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; } |