tests/test-fma2.h

   1 /* Test of fused multiply-add.
   2    Copyright (C) 2011 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
  16
  17 /* Written by Bruno Haible <bruno@clisp.org>, 2011.  */
  18
  19 /* Returns 2^e as a DOUBLE.  */
  20 #define POW2(e) \
  21   LDEXP (L_(1.0), e)
  22
  23 /* One could define XE_RANGE and YE_RANGE to 5 or 60, but this would slow down
  24    the test without the expectation of catching more bugs.  */
  25 #define XE_RANGE 0
  26 #define YE_RANGE 0
  27
  28 /* Define to 1 if you want to allow the behaviour of the 'double-double'
  29    implementation of 'long double' (seen on IRIX 6.5 and Linux/PowerPC).
  30    This floating-point type does not follow IEEE 754.  */
  31 #if MANT_BIT == LDBL_MANT_BIT && LDBL_MANT_BIT == 2 * DBL_MANT_BIT
  32 # define FORGIVE_DOUBLEDOUBLE_BUG 1
  33 #else
  34 # define FORGIVE_DOUBLEDOUBLE_BUG 0
  35 #endif
  36
  37 /* Subnormal numbers appear to not work as expected on IRIX 6.5.  */
  38 #ifdef __sgi
  39 # define MIN_SUBNORMAL_EXP (MIN_EXP - 1)
  40 #else
  41 # define MIN_SUBNORMAL_EXP (MIN_EXP - MANT_BIT)
  42 #endif
  43
  44 /* Check rounding behaviour.  */
  45
  46 static void
  47 test_function (DOUBLE (*my_fma) (DOUBLE, DOUBLE, DOUBLE))
  48 {
  49   /* Array mapping n to (-1)^n.  */
  50   static const DOUBLE pow_m1[] =
  51     {
  52       L_(1.0), - L_(1.0), L_(1.0), - L_(1.0),
  53       L_(1.0), - L_(1.0), L_(1.0), - L_(1.0)
  54     };
  55   volatile DOUBLE x;
  56   volatile DOUBLE y;
  57   volatile DOUBLE z;
  58   volatile DOUBLE result;
  59   volatile DOUBLE expected;
  60
  61   /* A product x * y that consists of two bits.  */
  62   {
  63     int xs;
  64     int xe;
  65     int ys;
  66     int ye;
  67     int ze;
  68     DOUBLE sign;
  69
  70     for (xs = 0; xs < 2; xs++)
  71       for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
  72         {
  73           x = pow_m1[xs] * POW2 (xe); /* (-1)^xs * 2^xe */
  74
  75           for (ys = 0; ys < 2; ys++)
  76             for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
  77               {
  78                 y = pow_m1[ys] * POW2 (ye); /* (-1)^ys * 2^ye */
  79
  80                 sign = pow_m1[xs + ys];
  81
  82                 /* Test addition (same signs).  */
  83                 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
  84                   {
  85                     z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
  86                     result = my_fma (x, y, z);
  87                     if (xe + ye >= ze + MANT_BIT)
  88                       expected = sign * POW2 (xe + ye);
  89                     else if (xe + ye > ze - MANT_BIT)
  90                       expected = sign * (POW2 (xe + ye) + POW2 (ze));
  91                     else
  92                       expected = z;
  93                     ASSERT (result == expected);
  94
  95                     ze++;
  96                     /* Shortcut some values of ze, to speed up the test.  */
  97                     if (ze == MIN_EXP + MANT_BIT)
  98                       ze = - 2 * MANT_BIT - 1;
  99                     else if (ze == 2 * MANT_BIT)
 100                       ze = MAX_EXP - MANT_BIT - 1;
 101                   }
 102
 103                 /* Test subtraction (opposite signs).  */
 104                 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
 105                   {
 106                     z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
 107                     result = my_fma (x, y, z);
 108                     if (xe + ye > ze + MANT_BIT)
 109                       expected = sign * POW2 (xe + ye);
 110                     else if (xe + ye >= ze)
 111                       expected = sign * (POW2 (xe + ye) - POW2 (ze));
 112                     else if (xe + ye > ze - 1 - MANT_BIT)
 113                       expected = - sign * (POW2 (ze) - POW2 (xe + ye));
 114                     else
 115                       expected = z;
 116                     ASSERT (result == expected);
 117
 118                     ze++;
 119                     /* Shortcut some values of ze, to speed up the test.  */
 120                     if (ze == MIN_EXP + MANT_BIT)
 121                       ze = - 2 * MANT_BIT - 1;
 122                     else if (ze == 2 * MANT_BIT)
 123                       ze = MAX_EXP - MANT_BIT - 1;
 124                   }
 125               }
 126         }
 127   }
 128   /* A product x * y that consists of three bits.  */
 129   {
 130     int i;
 131     int xs;
 132     int xe;
 133     int ys;
 134     int ye;
 135     int ze;
 136     DOUBLE sign;
 137
 138     for (i = 1; i <= MANT_BIT - 1; i++)
 139       for (xs = 0; xs < 2; xs++)
 140         for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
 141           {
 142             x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
 143               pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
 144
 145             for (ys = 0; ys < 2; ys++)
 146               for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
 147                 {
 148                   y = /* (-1)^ys * (2^ye + 2^(ye-i)) */
 149                     pow_m1[ys] * (POW2 (ye) + POW2 (ye - i));
 150
 151                   sign = pow_m1[xs + ys];
 152
 153                   /* The exact value of x * y is
 154                      (-1)^(xs+ys) * (2^(xe+ye) + 2^(xe+ye-i+1) + 2^(xe+ye-2*i)) */
 155
 156                   /* Test addition (same signs).  */
 157                   for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
 158                     {
 159                       z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
 160                       result = my_fma (x, y, z);
 161                       if (FORGIVE_DOUBLEDOUBLE_BUG)
 162                         if ((xe + ye > ze
 163                              && xe + ye < ze + MANT_BIT
 164                              && i == DBL_MANT_BIT)
 165                             || (xe + ye == ze + DBL_MANT_BIT && i == DBL_MANT_BIT + 1)
 166                             || (xe + ye == ze + MANT_BIT - 1 && i == 1))
 167                           goto skip1;
 168                       if (xe + ye > ze + MANT_BIT)
 169                         {
 170                           if (2 * i > MANT_BIT)
 171                             expected =
 172                               sign * (POW2 (xe + ye)
 173                                       + POW2 (xe + ye - i + 1));
 174                           else if (2 * i == MANT_BIT)
 175                             expected =
 176                               sign * (POW2 (xe + ye)
 177                                       + POW2 (xe + ye - i + 1)
 178                                       + POW2 (xe + ye - MANT_BIT + 1));
 179                           else
 180                             expected =
 181                               sign * (POW2 (xe + ye)
 182                                       + POW2 (xe + ye - i + 1)
 183                                       + POW2 (xe + ye - 2 * i));
 184                         }
 185                       else if (xe + ye == ze + MANT_BIT)
 186                         {
 187                           if (2 * i >= MANT_BIT)
 188                             expected =
 189                               sign * (POW2 (xe + ye)
 190                                       + POW2 (xe + ye - i + 1)
 191                                       + POW2 (xe + ye - MANT_BIT + 1));
 192                           else if (2 * i == MANT_BIT - 1)
 193                             /* round-to-even rounds up */
 194                             expected =
 195                               sign * (POW2 (xe + ye)
 196                                       + POW2 (xe + ye - i + 1)
 197                                       + POW2 (xe + ye - 2 * i + 1));
 198                           else
 199                             expected =
 200                               sign * (POW2 (xe + ye)
 201                                       + POW2 (xe + ye - i + 1)
 202                                       + POW2 (xe + ye - 2 * i));
 203                         }
 204                       else if (xe + ye > ze - MANT_BIT + 2 * i)
 205                         expected =
 206                           sign * (POW2 (ze)
 207                                   + POW2 (xe + ye)
 208                                   + POW2 (xe + ye - i + 1)
 209                                   + POW2 (xe + ye - 2 * i));
 210                       else if (xe + ye >= ze - MANT_BIT + i)
 211                         expected =
 212                           sign * (POW2 (ze)
 213                                   + POW2 (xe + ye)
 214                                   + POW2 (xe + ye - i + 1));
 215                       else if (xe + ye == ze - MANT_BIT + i - 1)
 216                         {
 217                           if (i == 1)
 218                             expected =
 219                               sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
 220                           else
 221                             expected =
 222                               sign * (POW2 (ze)
 223                                       + POW2 (xe + ye)
 224                                       + POW2 (ze - MANT_BIT + 1));
 225                         }
 226                       else if (xe + ye >= ze - MANT_BIT + 1)
 227                         expected = sign * (POW2 (ze) + POW2 (xe + ye));
 228                       else if (xe + ye == ze - MANT_BIT)
 229                         expected =
 230                           sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
 231                       else if (xe + ye == ze - MANT_BIT - 1)
 232                         {
 233                           if (i == 1)
 234                             expected =
 235                               sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
 236                           else
 237                             expected = z;
 238                         }
 239                       else
 240                         expected = z;
 241                       ASSERT (result == expected);
 242
 243                      skip1:
 244                       ze++;
 245                       /* Shortcut some values of ze, to speed up the test.  */
 246                       if (ze == MIN_EXP + MANT_BIT)
 247                         ze = - 2 * MANT_BIT - 1;
 248                       else if (ze == 2 * MANT_BIT)
 249                         ze = MAX_EXP - MANT_BIT - 1;
 250                     }
 251
 252                   /* Test subtraction (opposite signs).  */
 253                   if (i > 1)
 254                     for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
 255                       {
 256                         z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
 257                         result = my_fma (x, y, z);
 258                         if (FORGIVE_DOUBLEDOUBLE_BUG)
 259                           if ((xe + ye == ze && i == MANT_BIT - 1)
 260                               || (xe + ye > ze
 261                                   && xe + ye <= ze + DBL_MANT_BIT - 1
 262                                   && i == DBL_MANT_BIT + 1)
 263                               || (xe + ye >= ze + DBL_MANT_BIT - 1
 264                                   && xe + ye < ze + MANT_BIT
 265                                   && xe + ye == ze + i - 1)
 266                               || (xe + ye > ze + DBL_MANT_BIT
 267                                   && xe + ye < ze + MANT_BIT
 268                                   && i == DBL_MANT_BIT))
 269                             goto skip2;
 270                         if (xe + ye == ze)
 271                           {
 272                             /* maximal extinction */
 273                             expected =
 274                               sign * (POW2 (xe + ye - i + 1)
 275                                       + POW2 (xe + ye - 2 * i));
 276                           }
 277                         else if (xe + ye == ze - 1)
 278                           {
 279                             /* significant extinction */
 280                             if (2 * i > MANT_BIT)
 281                               expected =
 282                                 sign * (- POW2 (xe + ye)
 283                                         + POW2 (xe + ye - i + 1));
 284                             else
 285                               expected =
 286                                 sign * (- POW2 (xe + ye)
 287                                         + POW2 (xe + ye - i + 1)
 288                                         + POW2 (xe + ye - 2 * i));
 289                           }
 290                         else if (xe + ye > ze + MANT_BIT)
 291                           {
 292                             if (2 * i >= MANT_BIT)
 293                               expected =
 294                                 sign * (POW2 (xe + ye)
 295                                         + POW2 (xe + ye - i + 1));
 296                             else
 297                               expected =
 298                                 sign * (POW2 (xe + ye)
 299                                         + POW2 (xe + ye - i + 1)
 300                                         + POW2 (xe + ye - 2 * i));
 301                           }
 302                         else if (xe + ye == ze + MANT_BIT)
 303                           {
 304                             if (2 * i >= MANT_BIT)
 305                               expected =
 306                                 sign * (POW2 (xe + ye)
 307                                         + POW2 (xe + ye - i + 1));
 308                             else if (2 * i == MANT_BIT - 1)
 309                               /* round-to-even rounds down */
 310                               expected =
 311                                 sign * (POW2 (xe + ye)
 312                                         + POW2 (xe + ye - i + 1));
 313                             else
 314                               /* round-to-even rounds up */
 315                               expected =
 316                                 sign * (POW2 (xe + ye)
 317                                         + POW2 (xe + ye - i + 1)
 318                                         + POW2 (xe + ye - 2 * i));
 319                           }
 320                         else if (xe + ye >= ze - MANT_BIT + 2 * i)
 321                           expected =
 322                             sign * (- POW2 (ze)
 323                                     + POW2 (xe + ye)
 324                                     + POW2 (xe + ye - i + 1)
 325                                     + POW2 (xe + ye - 2 * i));
 326                         else if (xe + ye >= ze - MANT_BIT + i - 1)
 327                           expected =
 328                             sign * (- POW2 (ze)
 329                                     + POW2 (xe + ye)
 330                                     + POW2 (xe + ye - i + 1));
 331                         else if (xe + ye == ze - MANT_BIT + i - 2)
 332                           expected =
 333                             sign * (- POW2 (ze)
 334                                     + POW2 (xe + ye)
 335                                     + POW2 (ze - MANT_BIT));
 336                         else if (xe + ye >= ze - MANT_BIT)
 337                           expected =
 338                             sign * (- POW2 (ze)
 339                                     + POW2 (xe + ye));
 340                         else if (xe + ye == ze - MANT_BIT - 1)
 341                           expected =
 342                             sign * (- POW2 (ze)
 343                                     + POW2 (ze - MANT_BIT));
 344                         else
 345                           expected = z;
 346                         ASSERT (result == expected);
 347
 348                        skip2:
 349                         ze++;
 350                         /* Shortcut some values of ze, to speed up the test.  */
 351                         if (ze == MIN_EXP + MANT_BIT)
 352                           ze = - 2 * MANT_BIT - 1;
 353                         else if (ze == 2 * MANT_BIT)
 354                           ze = MAX_EXP - MANT_BIT - 1;
 355                       }
 356                 }
 357           }
 358   }
 359   /* TODO: Similar tests with
 360      x = (-1)^xs * (2^xe - 2^(xe-i)), y = (-1)^ys * (2^ye - 2^(ye-i))  */
 361   /* A product x * y that consists of one segment of bits (or, if you prefer,
 362      two bits, one with positive weight and one with negative weight).  */
 363   {
 364     int i;
 365     int xs;
 366     int xe;
 367     int ys;
 368     int ye;
 369     int ze;
 370     DOUBLE sign;
 371
 372     for (i = 1; i <= MANT_BIT - 1; i++)
 373       for (xs = 0; xs < 2; xs++)
 374         for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
 375           {
 376             x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
 377               pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
 378
 379             for (ys = 0; ys < 2; ys++)
 380               for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
 381                 {
 382                   y = /* (-1)^ys * (2^ye - 2^(ye-i)) */
 383                     pow_m1[ys] * (POW2 (ye) - POW2 (ye - i));
 384
 385                   sign = pow_m1[xs + ys];
 386
 387                   /* The exact value of x * y is
 388                      (-1)^(xs+ys) * (2^(xe+ye) - 2^(xe+ye-2*i)) */
 389
 390                   /* Test addition (same signs).  */
 391                   for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
 392                     {
 393                       z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
 394                       result = my_fma (x, y, z);
 395                       if (FORGIVE_DOUBLEDOUBLE_BUG)
 396                         if ((xe + ye == ze + MANT_BIT && i > DBL_MANT_BIT)
 397                             || (xe + ye < ze + MANT_BIT
 398                                 && xe + ye >= ze
 399                                 && i == DBL_MANT_BIT)
 400                             || (xe + ye < ze
 401                                 && xe + ye == ze - MANT_BIT + 2 * i))
 402                           goto skip3;
 403                       if (xe + ye > ze + MANT_BIT + 1)
 404                         {
 405                           if (2 * i > MANT_BIT)
 406                             expected = sign * POW2 (xe + ye);
 407                           else
 408                             expected =
 409                               sign * (POW2 (xe + ye)
 410                                       - POW2 (xe + ye - 2 * i));
 411                         }
 412                       else if (xe + ye == ze + MANT_BIT + 1)
 413                         {
 414                           if (2 * i >= MANT_BIT)
 415                             expected = sign * POW2 (xe + ye);
 416                           else
 417                             expected =
 418                               sign * (POW2 (xe + ye)
 419                                       - POW2 (xe + ye - 2 * i));
 420                         }
 421                       else if (xe + ye >= ze - MANT_BIT + 2 * i)
 422                         {
 423                           if (2 * i > MANT_BIT)
 424                             expected =
 425                               sign * (POW2 (xe + ye) + POW2 (ze));
 426                           else
 427                             expected =
 428                               sign * (POW2 (xe + ye)
 429                                       - POW2 (xe + ye - 2 * i)
 430                                       + POW2 (ze));
 431                         }
 432                       else if (xe + ye >= ze - MANT_BIT + 1)
 433                         expected =
 434                           sign * (POW2 (ze) + POW2 (xe + ye));
 435                       else
 436                         expected = z;
 437                       ASSERT (result == expected);
 438
 439                      skip3:
 440                       ze++;
 441                       /* Shortcut some values of ze, to speed up the test.  */
 442                       if (ze == MIN_EXP + MANT_BIT)
 443                         ze = - 2 * MANT_BIT - 1;
 444                       else if (ze == 2 * MANT_BIT)
 445                         ze = MAX_EXP - MANT_BIT - 1;
 446                     }
 447
 448                   /* Test subtraction (opposite signs).  */
 449                   if (i > 1)
 450                     for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
 451                       {
 452                         z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
 453                         result = my_fma (x, y, z);
 454                         if (FORGIVE_DOUBLEDOUBLE_BUG)
 455                           if (xe + ye > ze
 456                               && xe + ye < ze + DBL_MANT_BIT
 457                               && xe + ye == ze + 2 * i - LDBL_MANT_BIT)
 458                             goto skip4;
 459                         if (xe + ye == ze)
 460                           {
 461                             /* maximal extinction */
 462                             expected = sign * - POW2 (xe + ye - 2 * i);
 463                           }
 464                         else if (xe + ye > ze + MANT_BIT + 1)
 465                           {
 466                             if (2 * i > MANT_BIT + 1)
 467                               expected = sign * POW2 (xe + ye);
 468                             else if (2 * i == MANT_BIT + 1)
 469                               expected =
 470                                 sign * (POW2 (xe + ye)
 471                                         - POW2 (xe + ye - MANT_BIT));
 472                             else
 473                               expected =
 474                                 sign * (POW2 (xe + ye)
 475                                         - POW2 (xe + ye - 2 * i));
 476                           }
 477                         else if (xe + ye == ze + MANT_BIT + 1)
 478                           {
 479                             if (2 * i > MANT_BIT)
 480                               expected =
 481                                 sign * (POW2 (xe + ye)
 482                                         - POW2 (xe + ye - MANT_BIT));
 483                             else if (2 * i == MANT_BIT)
 484                               expected =
 485                                 sign * (POW2 (xe + ye)
 486                                         - POW2 (xe + ye - MANT_BIT + 1));
 487                             else
 488                               expected =
 489                                 sign * (POW2 (xe + ye)
 490                                         - POW2 (xe + ye - 2 * i));
 491                           }
 492                         else if (xe + ye == ze + MANT_BIT)
 493                           {
 494                             if (2 * i > MANT_BIT + 1)
 495                               expected =
 496                                 sign * (POW2 (xe + ye)
 497                                         - POW2 (xe + ye - MANT_BIT));
 498                             else if (2 * i == MANT_BIT + 1)
 499                               expected =
 500                                 sign * (POW2 (xe + ye)
 501                                         - POW2 (xe + ye - MANT_BIT + 1));
 502                             else
 503                               expected =
 504                                 sign * (POW2 (xe + ye)
 505                                         - POW2 (ze)
 506                                         - POW2 (xe + ye - 2 * i));
 507                           }
 508                         else if (xe + ye > ze - MANT_BIT + 2 * i)
 509                           {
 510                             if (2 * i > MANT_BIT)
 511                               expected =
 512                                 sign * (POW2 (xe + ye) - POW2 (ze));
 513                             else
 514                               expected =
 515                                 sign * (POW2 (xe + ye)
 516                                         - POW2 (ze)
 517                                         - POW2 (xe + ye - 2 * i));
 518                           }
 519                         else if (xe + ye == ze - MANT_BIT + 2 * i)
 520                           expected =
 521                             sign * (- POW2 (ze)
 522                                     + POW2 (xe + ye)
 523                                     - POW2 (xe + ye - 2 * i));
 524                         else if (xe + ye >= ze - MANT_BIT)
 525                           expected = sign * (- POW2 (ze) + POW2 (xe + ye));
 526                         else
 527                           expected = z;
 528                         ASSERT (result == expected);
 529
 530                        skip4:
 531                         ze++;
 532                         /* Shortcut some values of ze, to speed up the test.  */
 533                         if (ze == MIN_EXP + MANT_BIT)
 534                           ze = - 2 * MANT_BIT - 1;
 535                         else if (ze == 2 * MANT_BIT)
 536                           ze = MAX_EXP - MANT_BIT - 1;
 537                       }
 538                 }
 539           }
 540   }
 541   /* TODO: Tests with denormalized results.  */
 542   /* Tests with temporary overflow.  */
 543   {
 544     volatile DOUBLE x = POW2 (MAX_EXP - 1);
 545     volatile DOUBLE y = POW2 (MAX_EXP - 1);
 546     volatile DOUBLE z = - INFINITY;
 547     volatile DOUBLE result = my_fma (x, y, z);
 548     ASSERT (result == - INFINITY);
 549   }
 550   {
 551     volatile DOUBLE x = POW2 (MAX_EXP - 1); /* 2^(MAX_EXP-1) */
 552     volatile DOUBLE y = L_(2.0);            /* 2^1 */
 553     volatile DOUBLE z =               /* -(2^MAX_EXP - 2^(MAX_EXP-MANT_BIT)) */
 554       - LDEXP (POW2 (MAX_EXP - 1) - POW2 (MAX_EXP - MANT_BIT - 1), 1);
 555     volatile DOUBLE result = my_fma (x, y, z);
 556     if (!FORGIVE_DOUBLEDOUBLE_BUG)
 557       ASSERT (result == POW2 (MAX_EXP - MANT_BIT));
 558   }
 559   {
 560     volatile DOUBLE x = POW2 (MAX_EXP - 1);             /* 2^(MAX_EXP-1) */
 561     volatile DOUBLE y = L_(3.0);                        /* 3 */
 562     volatile DOUBLE z = - LDEXP (L_(5.0), MAX_EXP - 3); /* -5*2^(MAX_EXP-3) */
 563     volatile DOUBLE result = my_fma (x, y, z);
 564     ASSERT (result == LDEXP (L_(7.0), MAX_EXP - 3));
 565   }
 566 }