Start using fma in the libm implementation

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)  ; }