Tests for module 'fma'.

* modules/fma-tests: New file.
* tests/test-fma1.c: New file.
* tests/test-fma1.h: New file.
* tests/test-fma2.c: New file.
* tests/test-fma2.h: New file.

6 files changed, 971 insertions, 0 deletions

diff --git a/ChangeLog b/ChangeLog
index 3f5b617b37..d9af09c354 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,12 @@
 2011-11-05  Bruno Haible  <bruno@clisp.org>
 
+	Tests for module 'fma'.
+	* modules/fma-tests: New file.
+	* tests/test-fma1.c: New file.
+	* tests/test-fma1.h: New file.
+	* tests/test-fma2.c: New file.
+	* tests/test-fma2.h: New file.
+
 	New module 'fma'.
 	* lib/math.in.h (fma): New declaration.
 	* lib/fma.c: New file.
diff --git a/modules/fma-tests b/modules/fma-tests
new file mode 100644
index 0000000000..b1fbc7e346
--- /dev/null
+++ b/modules/fma-tests
@@ -0,0 +1,23 @@
+Files:
+tests/test-fma1.c
+tests/test-fma1.h
+tests/test-fma2.c
+tests/test-fma2.h
+tests/infinity.h
+tests/nan.h
+tests/signature.h
+tests/macros.h
+lib/float+.h
+
+Depends-on:
+float
+isnand-nolibm
+ldexp
+
+configure.ac:
+
+Makefile.am:
+TESTS += test-fma1 test-fma2
+check_PROGRAMS += test-fma1 test-fma2
+test_fma1_LDADD = $(LDADD) @FMA_LIBM@
+test_fma2_LDADD = $(LDADD) @FMA_LIBM@ @LDEXP_LIBM@
diff --git a/tests/test-fma1.c b/tests/test-fma1.c
new file mode 100644
index 0000000000..f3ef31e174
--- /dev/null
+++ b/tests/test-fma1.c
@@ -0,0 +1,47 @@
+/* Test of fma().
+   Copyright (C) 2011 Free Software Foundation, Inc.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation; either version 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+/* Written by Bruno Haible <bruno@clisp.org>, 2011.  */
+
+#include <config.h>
+
+#include <math.h>
+
+#include "signature.h"
+SIGNATURE_CHECK (fma, double, (double, double, double));
+
+#include "isnand-nolibm.h"
+#include "infinity.h"
+#include "nan.h"
+#include "macros.h"
+
+#undef INFINITY
+#undef NAN
+
+#define DOUBLE double
+#define ISNAN isnand
+#define INFINITY Infinityd ()
+#define NAN NaNd ()
+#define L_(literal) literal
+#include "test-fma1.h"
+
+int
+main ()
+{
+  test_function (fma);
+
+  return 0;
+}
diff --git a/tests/test-fma1.h b/tests/test-fma1.h
new file mode 100644
index 0000000000..3ea3c4a467
--- /dev/null
+++ b/tests/test-fma1.h
@@ -0,0 +1,213 @@
+/* Test of fused multiply-add.
+   Copyright (C) 2011 Free Software Foundation, Inc.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation; either version 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+/* Written by Bruno Haible <bruno@clisp.org>, 2011.  */
+
+static void
+test_function (DOUBLE (*my_fma) (DOUBLE, DOUBLE, DOUBLE))
+{
+  volatile DOUBLE x;
+  volatile DOUBLE y;
+  volatile DOUBLE z;
+  volatile DOUBLE result;
+  volatile DOUBLE expected;
+
+  /* Combinations with NaN.  */
+  /* "If x or y are NaN, a NaN shall be returned."  */
+  {
+    x = NAN;
+    y = NAN;
+    z = NAN;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = NAN;
+    y = NAN;
+    z = L_(1.0);
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = NAN;
+    y = L_(0.0);
+    z = NAN;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = NAN;
+    y = L_(0.0);
+    z = L_(1.0);
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = L_(0.0);
+    y = NAN;
+    z = NAN;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = L_(0.0);
+    y = NAN;
+    z = L_(1.0);
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  /* "If x*y is not 0*Inf nor Inf*0 and z is a NaN, a NaN shall be returned." */
+  {
+    x = L_(3.0);
+    y = - L_(2.0);
+    z = NAN;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  /* "If one of x and y is infinite, the other is zero, and z is a NaN, a NaN
+     shall be returned and a domain error may occur."  */
+  {
+    x = INFINITY;
+    y = L_(0.0);
+    z = NAN;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = L_(0.0);
+    y = INFINITY;
+    z = NAN;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+
+  /* Combinations with Infinity.  */
+  /* "If x multiplied by y is an exact infinity and z is also an infinity but
+     with the opposite sign, a domain error shall occur, and either a NaN
+     (if supported), or an implementation-defined value shall be returned."  */
+  {
+    x = INFINITY;
+    y = L_(3.0);
+    z = - INFINITY;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = INFINITY;
+    y = - L_(3.0);
+    z = INFINITY;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = L_(3.0);
+    y = INFINITY;
+    z = - INFINITY;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = - L_(3.0);
+    y = INFINITY;
+    z = INFINITY;
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  /* "If one of x and y is infinite, the other is zero, and z is not a NaN, a
+     domain error shall occur, and either a NaN (if supported), or an
+     implementation-defined value shall be returned."  */
+  {
+    x = INFINITY;
+    y = L_(0.0);
+    z = L_(5.0);
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  {
+    x = L_(0.0);
+    y = INFINITY;
+    z = L_(5.0);
+    result = my_fma (x, y, z);
+    ASSERT (ISNAN (result));
+  }
+  /* Infinite results.  */
+  {
+    x = - L_(2.0);
+    y = L_(3.0);
+    z = INFINITY;
+    result = my_fma (x, y, z);
+    expected = INFINITY;
+    ASSERT (result == expected);
+  }
+  {
+    x = INFINITY;
+    y = L_(3.0);
+    z = INFINITY;
+    result = my_fma (x, y, z);
+    expected = INFINITY;
+    ASSERT (result == expected);
+  }
+  {
+    x = INFINITY;
+    y = - L_(3.0);
+    z = - INFINITY;
+    result = my_fma (x, y, z);
+    expected = - INFINITY;
+    ASSERT (result == expected);
+  }
+  {
+    x = L_(3.0);
+    y = INFINITY;
+    z = INFINITY;
+    result = my_fma (x, y, z);
+    expected = INFINITY;
+    ASSERT (result == expected);
+  }
+  {
+    x = - L_(3.0);
+    y = INFINITY;
+    z = - INFINITY;
+    result = my_fma (x, y, z);
+    expected = - INFINITY;
+    ASSERT (result == expected);
+  }
+
+  /* Combinations with zero.  */
+  {
+    x = L_(0.0);
+    y = L_(3.0);
+    z = L_(11.0);
+    result = my_fma (x, y, z);
+    expected = L_(11.0);
+    ASSERT (result == expected);
+  }
+  {
+    x = L_(3.0);
+    y = L_(0.0);
+    z = L_(11.0);
+    result = my_fma (x, y, z);
+    expected = L_(11.0);
+    ASSERT (result == expected);
+  }
+  {
+    x = L_(3.0);
+    y = L_(4.0);
+    z = L_(0.0);
+    result = my_fma (x, y, z);
+    expected = L_(12.0);
+    ASSERT (result == expected);
+  }
+}
diff --git a/tests/test-fma2.c b/tests/test-fma2.c
new file mode 100644
index 0000000000..0fe91c7799
--- /dev/null
+++ b/tests/test-fma2.c
@@ -0,0 +1,47 @@
+/* Test of fma().
+   Copyright (C) 2011 Free Software Foundation, Inc.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation; either version 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+/* Written by Bruno Haible <bruno@clisp.org>, 2011.  */
+
+#include <config.h>
+
+#include <math.h>
+
+#include "float+.h"
+#include "infinity.h"
+#include "macros.h"
+
+#undef INFINITY
+
+#define DOUBLE double
+#define LDEXP ldexp
+const int MIN_EXP = DBL_MIN_EXP; /* for gdb */
+#define MIN_EXP DBL_MIN_EXP
+const int MAX_EXP = DBL_MAX_EXP; /* for gdb */
+#define MAX_EXP DBL_MAX_EXP
+const int MANT_BIT = DBL_MANT_BIT; /* for gdb */
+#define MANT_BIT DBL_MANT_BIT
+#define INFINITY Infinityd ()
+#define L_(literal) literal
+#include "test-fma2.h"
+
+int
+main ()
+{
+  test_function (fma);
+
+  return 0;
+}
diff --git a/tests/test-fma2.h b/tests/test-fma2.h
new file mode 100644
index 0000000000..b38bdbf53e
--- /dev/null
+++ b/tests/test-fma2.h
@@ -0,0 +1,634 @@
+/* Test of fused multiply-add.
+   Copyright (C) 2011 Free Software Foundation, Inc.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation; either version 3 of the License, or
+   (at your option) any later version.
+
+   This program 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 General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
+
+/* Written by Bruno Haible <bruno@clisp.org>, 2011.  */
+
+/* Returns 2^e as a DOUBLE.  */
+#define POW2(e) \
+  LDEXP (L_(1.0), e)
+
+/* One could define XE_RANGE and YE_RANGE to 5 or 60, but this would slow down
+   the test without the expectation of catching more bugs.  */
+#define XE_RANGE 0
+#define YE_RANGE 0
+
+/* Define to 1 if you want to allow the behaviour of the glibc 2.11 fma()
+   implementation.  glibc bug #1..#16 refer to the test cases in
+   <http://sourceware.org/bugzilla/show_bug.cgi?id=13304>.  */
+#if __GLIBC__ >= 2 && 0
+# define FORGIVE_GLIBC_BUG 1
+#else
+# define FORGIVE_GLIBC_BUG 0
+#endif
+
+/* Define to 1 if you want to allow the behaviour of the 'double-double'
+   implementation of 'long double' (seen on IRIX 6.5 and Linux/PowerPC).
+   This floating-point type does not follow IEEE 754.  */
+#if MANT_BIT == LDBL_MANT_BIT && LDBL_MANT_BIT == 2 * DBL_MANT_BIT
+# define FORGIVE_DOUBLEDOUBLE_BUG 1
+#else
+# define FORGIVE_DOUBLEDOUBLE_BUG 0
+#endif
+
+/* Subnormal numbers appear to not work as expected on IRIX 6.5.  */
+#ifdef __sgi
+# define MIN_SUBNORMAL_EXP (MIN_EXP - 1)
+#else
+# define MIN_SUBNORMAL_EXP (MIN_EXP - MANT_BIT)
+#endif
+
+/* Check rounding behaviour.  */
+
+static void
+test_function (DOUBLE (*my_fma) (DOUBLE, DOUBLE, DOUBLE))
+{
+  /* Array mapping n to (-1)^n.  */
+  static const DOUBLE pow_m1[] =
+    {
+      L_(1.0), - L_(1.0), L_(1.0), - L_(1.0),
+      L_(1.0), - L_(1.0), L_(1.0), - L_(1.0)
+    };
+  volatile DOUBLE x;
+  volatile DOUBLE y;
+  volatile DOUBLE z;
+  volatile DOUBLE result;
+  volatile DOUBLE expected;
+
+  /* A product x * y that consists of two bits.  */
+  {
+    int xs;
+    int xe;
+    int ys;
+    int ye;
+    int ze;
+    DOUBLE sign;
+
+    for (xs = 0; xs < 2; xs++)
+      for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
+        {
+          x = pow_m1[xs] * POW2 (xe); /* (-1)^xs * 2^xe */
+
+          for (ys = 0; ys < 2; ys++)
+            for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
+              {
+                y = pow_m1[ys] * POW2 (ye); /* (-1)^ys * 2^ye */
+
+                sign = pow_m1[xs + ys];
+
+                /* Test addition (same signs).  */
+                for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
+                  {
+                    z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
+                    result = my_fma (x, y, z);
+                    if (xe + ye >= ze + MANT_BIT)
+                      expected = sign * POW2 (xe + ye);
+                    else if (xe + ye > ze - MANT_BIT)
+                      expected = sign * (POW2 (xe + ye) + POW2 (ze));
+                    else
+                      expected = z;
+                    ASSERT (result == expected);
+
+                    ze++;
+                    /* Shortcut some values of ze, to speed up the test.  */
+                    if (ze == MIN_EXP + MANT_BIT)
+                      ze = - 2 * MANT_BIT - 1;
+                    else if (ze == 2 * MANT_BIT)
+                      ze = MAX_EXP - MANT_BIT - 1;
+                  }
+
+                /* Test subtraction (opposite signs).  */
+                for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
+                  {
+                    z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
+                    result = my_fma (x, y, z);
+                    if (xe + ye > ze + MANT_BIT)
+                      expected = sign * POW2 (xe + ye);
+                    else if (xe + ye >= ze)
+                      expected = sign * (POW2 (xe + ye) - POW2 (ze));
+                    else if (xe + ye > ze - 1 - MANT_BIT)
+                      expected = - sign * (POW2 (ze) - POW2 (xe + ye));
+                    else
+                      expected = z;
+                    ASSERT (result == expected);
+
+                    ze++;
+                    /* Shortcut some values of ze, to speed up the test.  */
+                    if (ze == MIN_EXP + MANT_BIT)
+                      ze = - 2 * MANT_BIT - 1;
+                    else if (ze == 2 * MANT_BIT)
+                      ze = MAX_EXP - MANT_BIT - 1;
+                  }
+              }
+        }
+  }
+  /* A product x * y that consists of three bits.  */
+  {
+    int i;
+    int xs;
+    int xe;
+    int ys;
+    int ye;
+    int ze;
+    DOUBLE sign;
+
+    for (i = 1; i <= MANT_BIT - 1; i++)
+      for (xs = 0; xs < 2; xs++)
+        for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
+          {
+            x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
+              pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
+
+            for (ys = 0; ys < 2; ys++)
+              for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
+                {
+                  y = /* (-1)^ys * (2^ye + 2^(ye-i)) */
+                    pow_m1[ys] * (POW2 (ye) + POW2 (ye - i));
+
+                  sign = pow_m1[xs + ys];
+
+                  /* The exact value of x * y is
+                     (-1)^(xs+ys) * (2^(xe+ye) + 2^(xe+ye-i+1) + 2^(xe+ye-2*i)) */
+
+                  /* Test addition (same signs).  */
+                  for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
+                    {
+                      z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
+                      result = my_fma (x, y, z);
+                      if (FORGIVE_DOUBLEDOUBLE_BUG)
+                        if ((xe + ye > ze
+                             && xe + ye < ze + MANT_BIT
+                             && i == DBL_MANT_BIT)
+                            || (xe + ye == ze + DBL_MANT_BIT && i == DBL_MANT_BIT + 1)
+                            || (xe + ye == ze + MANT_BIT - 1 && i == 1))
+                          goto skip1;
+                      if (xe + ye > ze + MANT_BIT)
+                        {
+                          if (2 * i > MANT_BIT)
+                            expected =
+                              sign * (POW2 (xe + ye)
+                                      + POW2 (xe + ye - i + 1));
+                          else if (2 * i == MANT_BIT)
+                            {
+                              if (FORGIVE_GLIBC_BUG) /* glibc bug #7 */
+                                goto skip1;
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        + POW2 (xe + ye - i + 1)
+                                        + POW2 (xe + ye - MANT_BIT + 1));
+                            }
+                          else
+                            expected =
+                              sign * (POW2 (xe + ye)
+                                      + POW2 (xe + ye - i + 1)
+                                      + POW2 (xe + ye - 2 * i));
+                        }
+                      else if (xe + ye == ze + MANT_BIT)
+                        {
+                          if (2 * i >= MANT_BIT)
+                            {
+                              if (FORGIVE_GLIBC_BUG) /* glibc bug #8 */
+                                goto skip1;
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        + POW2 (xe + ye - i + 1)
+                                        + POW2 (xe + ye - MANT_BIT + 1));
+                            }
+                          else if (2 * i == MANT_BIT - 1)
+                            /* round-to-even rounds up */
+                            expected =
+                              sign * (POW2 (xe + ye)
+                                      + POW2 (xe + ye - i + 1)
+                                      + POW2 (xe + ye - 2 * i + 1));
+                          else
+                            expected =
+                              sign * (POW2 (xe + ye)
+                                      + POW2 (xe + ye - i + 1)
+                                      + POW2 (xe + ye - 2 * i));
+                        }
+                      else if (xe + ye > ze - MANT_BIT + 2 * i)
+                        {
+                          if (2 * i == MANT_BIT)
+                            if (FORGIVE_GLIBC_BUG) /* glibc bug #9 */
+                              goto skip1;
+                          expected =
+                            sign * (POW2 (ze)
+                                    + POW2 (xe + ye)
+                                    + POW2 (xe + ye - i + 1)
+                                    + POW2 (xe + ye - 2 * i));
+                        }
+                      else if (xe + ye >= ze - MANT_BIT + i)
+                        expected =
+                          sign * (POW2 (ze)
+                                  + POW2 (xe + ye)
+                                  + POW2 (xe + ye - i + 1));
+                      else if (xe + ye == ze - MANT_BIT + i - 1)
+                        {
+                          if (2 * i >= MANT_BIT)
+                            if (FORGIVE_GLIBC_BUG) /* glibc bug #3, #10 */
+                              goto skip1;
+                          if (i == 1)
+                            expected =
+                              sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
+                          else
+                            expected =
+                              sign * (POW2 (ze)
+                                      + POW2 (xe + ye)
+                                      + POW2 (ze - MANT_BIT + 1));
+                        }
+                      else if (xe + ye >= ze - MANT_BIT + 1)
+                        expected = sign * (POW2 (ze) + POW2 (xe + ye));
+                      else if (xe + ye == ze - MANT_BIT)
+                        expected =
+                          sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
+                      else if (xe + ye == ze - MANT_BIT - 1)
+                        {
+                          if (i == 1)
+                            {
+                              if (FORGIVE_GLIBC_BUG) /* glibc bug #1 */
+                                goto skip1;
+                              expected =
+                                sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
+                            }
+                          else
+                            expected = z;
+                        }
+                      else
+                        expected = z;
+                      ASSERT (result == expected);
+
+                     skip1:
+                      ze++;
+                      /* Shortcut some values of ze, to speed up the test.  */
+                      if (ze == MIN_EXP + MANT_BIT)
+                        ze = - 2 * MANT_BIT - 1;
+                      else if (ze == 2 * MANT_BIT)
+                        ze = MAX_EXP - MANT_BIT - 1;
+                    }
+
+                  /* Test subtraction (opposite signs).  */
+                  if (i > 1)
+                    for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
+                      {
+                        z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
+                        result = my_fma (x, y, z);
+                        if (FORGIVE_DOUBLEDOUBLE_BUG)
+                          if ((xe + ye == ze && i == MANT_BIT - 1)
+                              || (xe + ye > ze
+                                  && xe + ye <= ze + DBL_MANT_BIT - 1
+                                  && i == DBL_MANT_BIT + 1)
+                              || (xe + ye >= ze + DBL_MANT_BIT - 1
+                                  && xe + ye < ze + MANT_BIT
+                                  && xe + ye == ze + i - 1)
+                              || (xe + ye > ze + DBL_MANT_BIT
+                                  && xe + ye < ze + MANT_BIT
+                                  && i == DBL_MANT_BIT))
+                            goto skip2;
+                        if (xe + ye == ze)
+                          {
+                            /* maximal extinction */
+                            if (2 * i >= MANT_BIT)
+                              if (FORGIVE_GLIBC_BUG) /* glibc bug #12 */
+                                goto skip2;
+                            expected =
+                              sign * (POW2 (xe + ye - i + 1)
+                                      + POW2 (xe + ye - 2 * i));
+                          }
+                        else if (xe + ye == ze - 1)
+                          {
+                            /* significant extinction */
+                            if (2 * i > MANT_BIT)
+                              expected =
+                                sign * (- POW2 (xe + ye)
+                                        + POW2 (xe + ye - i + 1));
+                            else
+                              {
+                                if (2 * i == MANT_BIT)
+                                  if (FORGIVE_GLIBC_BUG) /* glibc bug #13 */
+                                    goto skip2;
+                                expected =
+                                  sign * (- POW2 (xe + ye)
+                                          + POW2 (xe + ye - i + 1)
+                                          + POW2 (xe + ye - 2 * i));
+                              }
+                          }
+                        else if (xe + ye > ze + MANT_BIT)
+                          {
+                            if (2 * i >= MANT_BIT)
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        + POW2 (xe + ye - i + 1));
+                            else
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        + POW2 (xe + ye - i + 1)
+                                        + POW2 (xe + ye - 2 * i));
+                          }
+                        else if (xe + ye == ze + MANT_BIT)
+                          {
+                            if (2 * i >= MANT_BIT)
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        + POW2 (xe + ye - i + 1));
+                            else if (2 * i == MANT_BIT - 1)
+                              /* round-to-even rounds down */
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        + POW2 (xe + ye - i + 1));
+                            else
+                              /* round-to-even rounds up */
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        + POW2 (xe + ye - i + 1)
+                                        + POW2 (xe + ye - 2 * i));
+                          }
+                        else if (xe + ye >= ze - MANT_BIT + 2 * i)
+                          {
+                            if (2 * i == MANT_BIT)
+                              if (FORGIVE_GLIBC_BUG) /* glibc bug #11 */
+                                goto skip2;
+                            expected =
+                              sign * (- POW2 (ze)
+                                      + POW2 (xe + ye)
+                                      + POW2 (xe + ye - i + 1)
+                                      + POW2 (xe + ye - 2 * i));
+                          }
+                        else if (xe + ye >= ze - MANT_BIT + i - 1)
+                          expected =
+                            sign * (- POW2 (ze)
+                                    + POW2 (xe + ye)
+                                    + POW2 (xe + ye - i + 1));
+                        else if (xe + ye == ze - MANT_BIT + i - 2)
+                          {
+                            if (2 * i >= MANT_BIT)
+                              if (FORGIVE_GLIBC_BUG) /* glibc bug #4, #14 */
+                                goto skip2;
+                            expected =
+                              sign * (- POW2 (ze)
+                                      + POW2 (xe + ye)
+                                      + POW2 (ze - MANT_BIT));
+                          }
+                        else if (xe + ye >= ze - MANT_BIT)
+                          expected =
+                            sign * (- POW2 (ze)
+                                    + POW2 (xe + ye));
+                        else if (xe + ye == ze - MANT_BIT - 1)
+                          {
+                            if (FORGIVE_GLIBC_BUG) /* glibc bug #2 */
+                              goto skip2;
+                            expected =
+                              sign * (- POW2 (ze)
+                                      + POW2 (ze - MANT_BIT));
+                          }
+                        else
+                          expected = z;
+                        ASSERT (result == expected);
+
+                       skip2:
+                        ze++;
+                        /* Shortcut some values of ze, to speed up the test.  */
+                        if (ze == MIN_EXP + MANT_BIT)
+                          ze = - 2 * MANT_BIT - 1;
+                        else if (ze == 2 * MANT_BIT)
+                          ze = MAX_EXP - MANT_BIT - 1;
+                      }
+                }
+          }
+  }
+  /* TODO: Similar tests with
+     x = (-1)^xs * (2^xe - 2^(xe-i)), y = (-1)^ys * (2^ye - 2^(ye-i))  */
+  /* A product x * y that consists of one segment of bits (or, if you prefer,
+     two bits, one with positive weight and one with negative weight).  */
+  {
+    int i;
+    int xs;
+    int xe;
+    int ys;
+    int ye;
+    int ze;
+    DOUBLE sign;
+
+    for (i = 1; i <= MANT_BIT - 1; i++)
+      for (xs = 0; xs < 2; xs++)
+        for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
+          {
+            x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
+              pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
+
+            for (ys = 0; ys < 2; ys++)
+              for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
+                {
+                  y = /* (-1)^ys * (2^ye - 2^(ye-i)) */
+                    pow_m1[ys] * (POW2 (ye) - POW2 (ye - i));
+
+                  sign = pow_m1[xs + ys];
+
+                  /* The exact value of x * y is
+                     (-1)^(xs+ys) * (2^(xe+ye) - 2^(xe+ye-2*i)) */
+
+                  /* Test addition (same signs).  */
+                  for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
+                    {
+                      z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
+                      result = my_fma (x, y, z);
+                      if (FORGIVE_DOUBLEDOUBLE_BUG)
+                        if ((xe + ye == ze + MANT_BIT && i > DBL_MANT_BIT)
+                            || (xe + ye < ze + MANT_BIT
+                                && xe + ye >= ze
+                                && i == DBL_MANT_BIT)
+                            || (xe + ye < ze
+                                && xe + ye == ze - MANT_BIT + 2 * i))
+                          goto skip3;
+                      if (xe + ye > ze + MANT_BIT + 1)
+                        {
+                          if (2 * i > MANT_BIT)
+                            expected = sign * POW2 (xe + ye);
+                          else
+                            expected =
+                              sign * (POW2 (xe + ye)
+                                      - POW2 (xe + ye - 2 * i));
+                        }
+                      else if (xe + ye == ze + MANT_BIT + 1)
+                        {
+                          if (2 * i >= MANT_BIT)
+                            expected = sign * POW2 (xe + ye);
+                          else
+                            expected =
+                              sign * (POW2 (xe + ye)
+                                      - POW2 (xe + ye - 2 * i));
+                        }
+                      else if (xe + ye >= ze - MANT_BIT + 2 * i)
+                        {
+                          if (2 * i > MANT_BIT)
+                            expected =
+                              sign * (POW2 (xe + ye) + POW2 (ze));
+                          else
+                            expected =
+                              sign * (POW2 (xe + ye)
+                                      - POW2 (xe + ye - 2 * i)
+                                      + POW2 (ze));
+                        }
+                      else if (xe + ye >= ze - MANT_BIT + 1)
+                        expected =
+                          sign * (POW2 (ze) + POW2 (xe + ye));
+                      else
+                        expected = z;
+                      ASSERT (result == expected);
+
+                     skip3:
+                      ze++;
+                      /* Shortcut some values of ze, to speed up the test.  */
+                      if (ze == MIN_EXP + MANT_BIT)
+                        ze = - 2 * MANT_BIT - 1;
+                      else if (ze == 2 * MANT_BIT)
+                        ze = MAX_EXP - MANT_BIT - 1;
+                    }
+
+                  /* Test subtraction (opposite signs).  */
+                  if (i > 1)
+                    for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
+                      {
+                        z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
+                        result = my_fma (x, y, z);
+                        if (FORGIVE_DOUBLEDOUBLE_BUG)
+                          if (xe + ye > ze
+                              && xe + ye < ze + DBL_MANT_BIT
+                              && xe + ye == ze + 2 * i - LDBL_MANT_BIT)
+                            goto skip4;
+                        if (xe + ye == ze)
+                          {
+                            /* maximal extinction */
+                            if (2 * i > MANT_BIT)
+                              if (FORGIVE_GLIBC_BUG) /* glibc bug #16 */
+                                goto skip4;
+                            expected = sign * - POW2 (xe + ye - 2 * i);
+                          }
+                        else if (xe + ye > ze + MANT_BIT + 1)
+                          {
+                            if (2 * i > MANT_BIT + 1)
+                              expected = sign * POW2 (xe + ye);
+                            else if (2 * i == MANT_BIT + 1)
+                              {
+                                if (FORGIVE_GLIBC_BUG) /* glibc bug #6 */
+                                  goto skip4;
+                                expected =
+                                  sign * (POW2 (xe + ye)
+                                          - POW2 (xe + ye - MANT_BIT));
+                              }
+                            else
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        - POW2 (xe + ye - 2 * i));
+                          }
+                        else if (xe + ye == ze + MANT_BIT + 1)
+                          {
+                            if (2 * i > MANT_BIT)
+                              {
+                                if (FORGIVE_GLIBC_BUG) /* glibc bug #15 */
+                                  goto skip4;
+                                expected =
+                                  sign * (POW2 (xe + ye)
+                                          - POW2 (xe + ye - MANT_BIT));
+                              }
+                            else if (2 * i == MANT_BIT)
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        - POW2 (xe + ye - MANT_BIT + 1));
+                            else
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        - POW2 (xe + ye - 2 * i));
+                          }
+                        else if (xe + ye == ze + MANT_BIT)
+                          {
+                            if (2 * i > MANT_BIT + 1)
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        - POW2 (xe + ye - MANT_BIT));
+                            else if (2 * i == MANT_BIT + 1)
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        - POW2 (xe + ye - MANT_BIT + 1));
+                            else
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        - POW2 (ze)
+                                        - POW2 (xe + ye - 2 * i));
+                          }
+                        else if (xe + ye > ze - MANT_BIT + 2 * i)
+                          {
+                            if (2 * i > MANT_BIT)
+                              expected =
+                                sign * (POW2 (xe + ye) - POW2 (ze));
+                            else
+                              expected =
+                                sign * (POW2 (xe + ye)
+                                        - POW2 (ze)
+                                        - POW2 (xe + ye - 2 * i));
+                          }
+                        else if (xe + ye == ze - MANT_BIT + 2 * i)
+                          expected =
+                            sign * (- POW2 (ze)
+                                    + POW2 (xe + ye)
+                                    - POW2 (xe + ye - 2 * i));
+                        else if (xe + ye > ze - MANT_BIT)
+                          expected = sign * (- POW2 (ze) + POW2 (xe + ye));
+                        else if (xe + ye == ze - MANT_BIT)
+                          {
+                            if (FORGIVE_GLIBC_BUG) /* glibc bug #5 */
+                              goto skip4;
+                            expected = sign * (- POW2 (ze) + POW2 (xe + ye));
+                          }
+                        else
+                          expected = z;
+                        ASSERT (result == expected);
+
+                       skip4:
+                        ze++;
+                        /* Shortcut some values of ze, to speed up the test.  */
+                        if (ze == MIN_EXP + MANT_BIT)
+                          ze = - 2 * MANT_BIT - 1;
+                        else if (ze == 2 * MANT_BIT)
+                          ze = MAX_EXP - MANT_BIT - 1;
+                      }
+                }
+          }
+  }
+  /* TODO: Tests with denormalized results.  */
+  /* Tests with temporary overflow.  */
+  {
+    volatile DOUBLE x = POW2 (MAX_EXP - 1);
+    volatile DOUBLE y = POW2 (MAX_EXP - 1);
+    volatile DOUBLE z = - INFINITY;
+    volatile DOUBLE result = my_fma (x, y, z);
+    ASSERT (result == - INFINITY);
+  }
+  {
+    volatile DOUBLE x = POW2 (MAX_EXP - 1); /* 2^(MAX_EXP-1) */
+    volatile DOUBLE y = L_(2.0);            /* 2^1 */
+    volatile DOUBLE z =               /* -(2^MAX_EXP - 2^(MAX_EXP-MANT_BIT)) */
+      - LDEXP (POW2 (MAX_EXP - 1) - POW2 (MAX_EXP - MANT_BIT - 1), 1);
+    volatile DOUBLE result = my_fma (x, y, z);
+    if (!FORGIVE_DOUBLEDOUBLE_BUG)
+      ASSERT (result == POW2 (MAX_EXP - MANT_BIT));
+  }
+  {
+    volatile DOUBLE x = POW2 (MAX_EXP - 1);             /* 2^(MAX_EXP-1) */
+    volatile DOUBLE y = L_(3.0);                        /* 3 */
+    volatile DOUBLE z = - LDEXP (L_(5.0), MAX_EXP - 3); /* -5*2^(MAX_EXP-3) */
+    volatile DOUBLE result = my_fma (x, y, z);
+    ASSERT (result == LDEXP (L_(7.0), MAX_EXP - 3));
+  }
+}