Mercurial > gnulib
changeset 16039:ae7cb951d5b4
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.
| author | Bruno Haible <bruno@clisp.org> |
|---|---|
| date | Sun, 06 Nov 2011 02:28:32 +0100 |
| parents | 8f14d00d6f24 |
| children | b8acd8099b25 |
| files | ChangeLog modules/fma-tests tests/test-fma1.c tests/test-fma1.h tests/test-fma2.c tests/test-fma2.h |
| diffstat | 6 files changed, 971 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog Mon Oct 17 23:48:01 2011 +0200 +++ b/ChangeLog Sun Nov 06 02:28:32 2011 +0100 @@ -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.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/modules/fma-tests Sun Nov 06 02:28:32 2011 +0100 @@ -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@
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/test-fma1.c Sun Nov 06 02:28:32 2011 +0100 @@ -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; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/test-fma1.h Sun Nov 06 02:28:32 2011 +0100 @@ -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); + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/test-fma2.c Sun Nov 06 02:28:32 2011 +0100 @@ -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; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/test-fma2.h Sun Nov 06 02:28:32 2011 +0100 @@ -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)); + } +}