This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[PATCH][BZ #3268] Add fma single/double implementations to soft-fp


The soft-fp implementation of IBM extended long double exposed a general
lack to correct fma implemetations for platforms that don't implement
fma in hardware (soft-fp plus hard-fp that doesn't implement fma).

This patch provides a basic soft-fp implementation for single and double
fma.  I assume that __fmadf4/__fmasf4 are acceptable names for these
functions. These implementation needs to be in libc.so to access the
soft-float exeception and rounding modes. Then the s_fma.c/s_fmaf.c (in
libm.so) need to be overridden for nofpu to call __fmadf4/__fmasf4 directly.

This patch also updates libm-test.inc to actually test for the full
percision of fma.

2006-10-05  Steven Munroe  <sjmunroe@us.ibm.com>

	[BZ #3268]
	* soft-fp/Makefile (gcc-single-routines): Add fmasf4.
	(gcc-double-routines): Add fmadf4.
	* soft-fp/fmadf4.c: New file.
	* soft-fp/fmasf4.c: New file.

	[BZ #3268]
	* math/libm-test.inc (fma_test): New tests.

diff -urN libc25-cvstip-20061005/soft-fp/Makefile libc24/soft-fp/Makefile
--- libc25-cvstip-20061005/soft-fp/Makefile	2006-01-06 04:47:45.000000000 -0600
+++ libc24/soft-fp/Makefile	2006-10-10 09:54:48.000000000 -0500
@@ -24,12 +24,13 @@
 
 gcc-single-routines := negsf2 addsf3 subsf3 mulsf3 divsf3 eqsf2 \
 	lesf2 gesf2 unordsf2 fixsfsi fixunssfsi floatsisf fixsfdi        \
-	fixunssfdi floatdisf sqrtsf2 floatunsisf floatundisf
+	fixunssfdi floatdisf sqrtsf2 floatunsisf floatundisf \
+	fmasf4
 
 gcc-double-routines := negdf2 adddf3 subdf3 muldf3 divdf3 eqdf2 \
 	ledf2 gedf2 unorddf2 fixdfsi fixunsdfsi floatsidf fixdfdi        \
 	fixunsdfdi floatdidf extendsfdf2 truncdfsf2 sqrtdf2 floatunsidf \
-	floatundidf
+	floatundidf fmadf4
 
 gcc-quad-routines := negtf2 addtf3 subtf3 multf3 divtf3 eqtf2 \
 	letf2 getf2 unordtf2 fixtfsi fixunstfsi floatsitf fixtfdi      \
diff -urN libc25-cvstip-20061005/soft-fp/fmadf4.c libc24/soft-fp/fmadf4.c
--- libc25-cvstip-20061005/soft-fp/fmadf4.c	Wed Dec 31 18:00:00 1969
+++ libc24/soft-fp/fmadf4.c	Tue Oct 10 09:54:48 2006
@@ -0,0 +1,97 @@
+/* soft-fp x * y + z as ternary operation.
+   Copyright (C) 2006 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Steven Munroe <sjmunroe@us.ibm.com>, 2006.
+
+   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.
+
+   In addition to the permissions in the GNU Lesser General Public
+   License, the Free Software Foundation gives you unlimited
+   permission to link the compiled version of this file into
+   combinations with other programs, and to distribute those
+   combinations without any restriction coming from the use of this
+   file.  (The Lesser General Public License restrictions do apply in
+   other respects; for example, they cover modification of the file,
+   and distribution when not linked into a combine executable.)
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include "soft-fp.h"
+#include "double.h"
+#include "quad.h"
+
+/* Compute floating point multiply-add with higher (quad) precision.  */
+DFtype
+__fmadf4 (DFtype a, DFtype b, DFtype c)
+{
+    FP_DECL_EX;
+    FP_DECL_D(A);
+    FP_DECL_D(B);
+    FP_DECL_D(C);
+    FP_DECL_Q(X);
+    FP_DECL_Q(Y);
+    FP_DECL_Q(Z);
+    FP_DECL_Q(U);
+    FP_DECL_Q(V);
+    FP_DECL_D(R);
+    double r;
+    long double u, x, y, z;
+  
+    FP_INIT_ROUNDMODE;
+    FP_UNPACK_RAW_D (A, a);
+    FP_UNPACK_RAW_D (B, b);
+    FP_UNPACK_RAW_D (C, c);
+
+    /* Extend double to quad.  */
+#if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
+    FP_EXTEND(Q,D,4,2,X,A);
+    FP_EXTEND(Q,D,4,2,Y,B);
+    FP_EXTEND(Q,D,4,2,Z,C);
+#else
+    FP_EXTEND(Q,D,2,1,X,A);
+    FP_EXTEND(Q,D,2,1,Y,B);
+    FP_EXTEND(Q,D,2,1,Z,C);
+#endif
+    FP_PACK_RAW_Q(x,X);
+    FP_PACK_RAW_Q(y,Y);
+    FP_PACK_RAW_Q(z,Z);
+    FP_HANDLE_EXCEPTIONS;
+
+    /* Multiply.
+       Rounding is not an issue as we keep the full 106 bit product.  */
+    FP_UNPACK_Q(X,x);
+    FP_UNPACK_Q(Y,y);
+    FP_MUL_Q(U,X,Y);
+    FP_PACK_Q(u,U);
+    FP_HANDLE_EXCEPTIONS;
+
+    /* Add without rounding.  */
+    FP_UNPACK_SEMIRAW_Q(U,u);
+    FP_UNPACK_SEMIRAW_Q(Z,z);
+    FP_ADD_Q(V,U,Z);
+
+    /* Truncate quad to double and round.  */
+#if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
+    V_f[3] &= 0x0007ffff;  /* SEMIRAW_TO_SEMIRAW(4) */
+    FP_TRUNC(D,Q,2,4,R,V);
+#else
+    V_f1 &= 0x0007ffffffffffffL;  /* SEMIRAW_TO_SEMIRAW(2) */
+    FP_TRUNC(D,Q,1,2,R,V);
+#endif
+    FP_PACK_SEMIRAW_D(r,R);
+    FP_HANDLE_EXCEPTIONS;
+
+    return r;
+}
+
diff -urN libc25-cvstip-20061005/soft-fp/fmasf4.c libc24/soft-fp/fmasf4.c
--- libc25-cvstip-20061005/soft-fp/fmasf4.c	Wed Dec 31 18:00:00 1969
+++ libc24/soft-fp/fmasf4.c	Tue Oct 10 09:54:48 2006
@@ -0,0 +1,96 @@
+/* soft-fp x * y + z as ternary operation.
+   Copyright (C) 2006 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Steven Munroe <sjmunroe@us.ibm.com>, 2006.
+
+   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.
+
+   In addition to the permissions in the GNU Lesser General Public
+   License, the Free Software Foundation gives you unlimited
+   permission to link the compiled version of this file into
+   combinations with other programs, and to distribute those
+   combinations without any restriction coming from the use of this
+   file.  (The Lesser General Public License restrictions do apply in
+   other respects; for example, they cover modification of the file,
+   and distribution when not linked into a combine executable.)
+
+   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, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include "soft-fp.h"
+#include "single.h"
+#include "double.h"
+
+/* Compute floating point multiply-add with higher (double) precision.  */
+SFtype
+__fmasf4 (SFtype a, SFtype b, SFtype c)
+{
+    FP_DECL_EX;
+    FP_DECL_S(A);
+    FP_DECL_S(B);
+    FP_DECL_S(C);
+    FP_DECL_D(X);
+    FP_DECL_D(Y);
+    FP_DECL_D(Z);
+    FP_DECL_D(U);
+    FP_DECL_D(V);
+    FP_DECL_S(R);
+    float r;
+    double u, x, y, z;
+
+    FP_INIT_ROUNDMODE;
+    FP_UNPACK_RAW_S (A, a);
+    FP_UNPACK_RAW_S (B, b);
+    FP_UNPACK_RAW_S (C, c);
+
+    /* Extend single to double.  */
+#if _FP_W_TYPE_SIZE < _FP_FRACBITS_D
+    FP_EXTEND(D,S,2,1,X,A);
+    FP_EXTEND(D,S,2,1,Y,B);
+    FP_EXTEND(D,S,2,1,Z,C);
+#else
+    FP_EXTEND(D,S,1,1,X,A);
+    FP_EXTEND(D,S,1,1,Y,B);
+    FP_EXTEND(D,S,1,1,Z,C);
+#endif
+    FP_PACK_RAW_D(x,X);
+    FP_PACK_RAW_D(y,Y);
+    FP_PACK_RAW_D(z,Z);
+    FP_HANDLE_EXCEPTIONS;
+
+    /* Multiply.
+       Rounding is not an issue as we keep the full 48 bit product.  */
+    FP_UNPACK_D(X,x);
+    FP_UNPACK_D(Y,y);
+    FP_MUL_D(U,X,Y);
+    FP_PACK_D(u,U);
+    FP_HANDLE_EXCEPTIONS;
+
+    /* Add without rounding.  */
+    FP_UNPACK_SEMIRAW_D(U,u);
+    FP_UNPACK_SEMIRAW_D(Z,z);
+    FP_ADD_D(V,U,Z);
+
+    /* Truncate double to single and round.  */
+#if FP_W_TYPE_SIZE < _FP_FRACBITS_D
+    V_f1 &= 0x007fffffL;  /* SEMIRAW_TO_SEMIRAW(2) */
+    FP_TRUNC(S,D,1,2,R,V);
+#else
+    V_f &= 0x007fffffffffffffL;  /* SEMIRAW_TO_SEMIRAW(1) */
+    FP_TRUNC(S,D,1,1,R,V);
+#endif
+    FP_PACK_SEMIRAW_S(r,R);
+    FP_HANDLE_EXCEPTIONS;
+
+    return r;
+}
diff -urN libc24-cvstip-20060920/math/libm-test.inc libc24/math/libm-test.inc
--- libc24-cvstip-20060920/math/libm-test.inc	2006-09-20 14:36:26.000000000 -0500
+++ libc24/math/libm-test.inc	2006-10-04 16:22:37.000000000 -0500
@@ -2766,6 +2766,12 @@
 
   TEST_fff_f (fma, 1.25L, 0.75L, 0.0625L, 1.0L);
 
+  TEST_fff_f (fma, 8388609.0L, 8388609.0L, -70368752566272.0L, 8388609.0L);
+
+#ifdef TEST_DOUBLE
+  TEST_fff_f (fma, 0x1.0000000000001p+52, 0x1.0000000000001p+52, -0x1.0000000000001p+104, 0x1.0000000000001p+52);
+#endif
+
   END (fma);
 }

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]