This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH][BZ #3268] Add fma single/double implementations to soft-fp
- From: Steven Munroe <munroesj at us dot ibm dot com>
- To: libc-alpha at sources dot redhat dot com, libc-ports at sources dot redhat dot com, David Edelsohn <dje at watson dot ibm dot com>, "Joseph S. Myers" <joseph at codesourcery dot com>
- Date: Thu, 12 Oct 2006 11:47:06 -0500
- Subject: [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);
}