/*
* 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
/* 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
/* 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. */
/* 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, */
/* 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. */
/* 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 */
/* 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 */
/* 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 */
/* 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 */
/* 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. */
/* 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. */
/* 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) ; }
projects / thirdparty / glibc.git / blobdiff