tests/test-fma2.h

   1 /* Test of fused multiply-add.
   2    Copyright (C) 2011-2013 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
  56   /* A product x * y that consists of two bits.  */
  57   {
  58     volatile DOUBLE x;
  59     volatile DOUBLE y;
  60     volatile DOUBLE z;
  61     volatile DOUBLE result;
  62     volatile DOUBLE expected;
  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     volatile DOUBLE x;
 131     volatile DOUBLE y;
 132     volatile DOUBLE z;
 133     volatile DOUBLE result;
 134     volatile DOUBLE expected;
 135     int i;
 136     int xs;
 137     int xe;
 138     int ys;
 139     int ye;
 140     int ze;
 141     DOUBLE sign;
 142
 143     for (i = 1; i <= MANT_BIT - 1; i++)
 144       for (xs = 0; xs < 2; xs++)
 145         for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
 146           {
 147             x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
 148               pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
 149
 150             for (ys = 0; ys < 2; ys++)
 151               for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
 152                 {
 153                   y = /* (-1)^ys * (2^ye + 2^(ye-i)) */
 154                     pow_m1[ys] * (POW2 (ye) + POW2 (ye - i));
 155
 156                   sign = pow_m1[xs + ys];
 157
 158                   /* The exact value of x * y is
 159                      (-1)^(xs+ys) * (2^(xe+ye) + 2^(xe+ye-i+1) + 2^(xe+ye-2*i)) */
 160
 161                   /* Test addition (same signs).  */
 162                   for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
 163                     {
 164                       z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
 165                       result = my_fma (x, y, z);
 166                       if (FORGIVE_DOUBLEDOUBLE_BUG)
 167                         if ((xe + ye > ze
 168                              && xe + ye < ze + MANT_BIT
 169                              && i == DBL_MANT_BIT)
 170                             || (xe + ye == ze + DBL_MANT_BIT && i == DBL_MANT_BIT + 1)
 171                             || (xe + ye == ze + MANT_BIT - 1 && i == 1))
 172                           goto skip1;
 173                       if (xe + ye > ze + MANT_BIT)
 174                         {
 175                           if (2 * i > MANT_BIT)
 176                             expected =
 177                               sign * (POW2 (xe + ye)
 178                                       + POW2 (xe + ye - i + 1));
 179                           else if (2 * i == MANT_BIT)
 180                             expected =
 181                               sign * (POW2 (xe + ye)
 182                                       + POW2 (xe + ye - i + 1)
 183                                       + POW2 (xe + ye - MANT_BIT + 1));
 184                           else
 185                             expected =
 186                               sign * (POW2 (xe + ye)
 187                                       + POW2 (xe + ye - i + 1)
 188                                       + POW2 (xe + ye - 2 * i));
 189                         }
 190                       else if (xe + ye == ze + MANT_BIT)
 191                         {
 192                           if (2 * i >= MANT_BIT)
 193                             expected =
 194                               sign * (POW2 (xe + ye)
 195                                       + POW2 (xe + ye - i + 1)
 196                                       + POW2 (xe + ye - MANT_BIT + 1));
 197                           else if (2 * i == MANT_BIT - 1)
 198                             /* round-to-even rounds up */
 199                             expected =
 200                               sign * (POW2 (xe + ye)
 201                                       + POW2 (xe + ye - i + 1)
 202                                       + POW2 (xe + ye - 2 * i + 1));
 203                           else
 204                             expected =
 205                               sign * (POW2 (xe + ye)
 206                                       + POW2 (xe + ye - i + 1)
 207                                       + POW2 (xe + ye - 2 * i));
 208                         }
 209                       else if (xe + ye > ze - MANT_BIT + 2 * i)
 210                         expected =
 211                           sign * (POW2 (ze)
 212                                   + POW2 (xe + ye)
 213                                   + POW2 (xe + ye - i + 1)
 214                                   + POW2 (xe + ye - 2 * i));
 215                       else if (xe + ye >= ze - MANT_BIT + i)
 216                         expected =
 217                           sign * (POW2 (ze)
 218                                   + POW2 (xe + ye)
 219                                   + POW2 (xe + ye - i + 1));
 220                       else if (xe + ye == ze - MANT_BIT + i - 1)
 221                         {
 222                           if (i == 1)
 223                             expected =
 224                               sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
 225                           else
 226                             expected =
 227                               sign * (POW2 (ze)
 228                                       + POW2 (xe + ye)
 229                                       + POW2 (ze - MANT_BIT + 1));
 230                         }
 231                       else if (xe + ye >= ze - MANT_BIT + 1)
 232                         expected = sign * (POW2 (ze) + POW2 (xe + ye));
 233                       else if (xe + ye == ze - MANT_BIT)
 234                         expected =
 235                           sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
 236                       else if (xe + ye == ze - MANT_BIT - 1)
 237                         {
 238                           if (i == 1)
 239                             expected =
 240                               sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
 241                           else
 242                             expected = z;
 243                         }
 244                       else
 245                         expected = z;
 246                       ASSERT (result == expected);
 247
 248                      skip1:
 249                       ze++;
 250                       /* Shortcut some values of ze, to speed up the test.  */
 251                       if (ze == MIN_EXP + MANT_BIT)
 252                         ze = - 2 * MANT_BIT - 1;
 253                       else if (ze == 2 * MANT_BIT)
 254                         ze = MAX_EXP - MANT_BIT - 1;
 255                     }
 256
 257                   /* Test subtraction (opposite signs).  */
 258                   if (i > 1)
 259                     for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
 260                       {
 261                         z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
 262                         result = my_fma (x, y, z);
 263                         if (FORGIVE_DOUBLEDOUBLE_BUG)
 264                           if ((xe + ye == ze && i == MANT_BIT - 1)
 265                               || (xe + ye > ze
 266                                   && xe + ye <= ze + DBL_MANT_BIT - 1
 267                                   && i == DBL_MANT_BIT + 1)
 268                               || (xe + ye >= ze + DBL_MANT_BIT - 1
 269                                   && xe + ye < ze + MANT_BIT
 270                                   && xe + ye == ze + i - 1)
 271                               || (xe + ye > ze + DBL_MANT_BIT
 272                                   && xe + ye < ze + MANT_BIT
 273                                   && i == DBL_MANT_BIT))
 274                             goto skip2;
 275                         if (xe + ye == ze)
 276                           {
 277                             /* maximal extinction */
 278                             expected =
 279                               sign * (POW2 (xe + ye - i + 1)
 280                                       + POW2 (xe + ye - 2 * i));
 281                           }
 282                         else if (xe + ye == ze - 1)
 283                           {
 284                             /* significant extinction */
 285                             if (2 * i > MANT_BIT)
 286                               expected =
 287                                 sign * (- POW2 (xe + ye)
 288                                         + POW2 (xe + ye - i + 1));
 289                             else
 290                               expected =
 291                                 sign * (- POW2 (xe + ye)
 292                                         + POW2 (xe + ye - i + 1)
 293                                         + POW2 (xe + ye - 2 * i));
 294                           }
 295                         else if (xe + ye > ze + MANT_BIT)
 296                           {
 297                             if (2 * i >= MANT_BIT)
 298                               expected =
 299                                 sign * (POW2 (xe + ye)
 300                                         + POW2 (xe + ye - i + 1));
 301                             else
 302                               expected =
 303                                 sign * (POW2 (xe + ye)
 304                                         + POW2 (xe + ye - i + 1)
 305                                         + POW2 (xe + ye - 2 * i));
 306                           }
 307                         else if (xe + ye == ze + MANT_BIT)
 308                           {
 309                             if (2 * i >= MANT_BIT)
 310                               expected =
 311                                 sign * (POW2 (xe + ye)
 312                                         + POW2 (xe + ye - i + 1));
 313                             else if (2 * i == MANT_BIT - 1)
 314                               /* round-to-even rounds down */
 315                               expected =
 316                                 sign * (POW2 (xe + ye)
 317                                         + POW2 (xe + ye - i + 1));
 318                             else
 319                               /* round-to-even rounds up */
 320                               expected =
 321                                 sign * (POW2 (xe + ye)
 322                                         + POW2 (xe + ye - i + 1)
 323                                         + POW2 (xe + ye - 2 * i));
 324                           }
 325                         else if (xe + ye >= ze - MANT_BIT + 2 * i)
 326                           expected =
 327                             sign * (- POW2 (ze)
 328                                     + POW2 (xe + ye)
 329                                     + POW2 (xe + ye - i + 1)
 330                                     + POW2 (xe + ye - 2 * i));
 331                         else if (xe + ye >= ze - MANT_BIT + i - 1)
 332                           expected =
 333                             sign * (- POW2 (ze)
 334                                     + POW2 (xe + ye)
 335                                     + POW2 (xe + ye - i + 1));
 336                         else if (xe + ye == ze - MANT_BIT + i - 2)
 337                           expected =
 338                             sign * (- POW2 (ze)
 339                                     + POW2 (xe + ye)
 340                                     + POW2 (ze - MANT_BIT));
 341                         else if (xe + ye >= ze - MANT_BIT)
 342                           expected =
 343                             sign * (- POW2 (ze)
 344                                     + POW2 (xe + ye));
 345                         else if (xe + ye == ze - MANT_BIT - 1)
 346                           expected =
 347                             sign * (- POW2 (ze)
 348                                     + POW2 (ze - MANT_BIT));
 349                         else
 350                           expected = z;
 351                         ASSERT (result == expected);
 352
 353                        skip2:
 354                         ze++;
 355                         /* Shortcut some values of ze, to speed up the test.  */
 356                         if (ze == MIN_EXP + MANT_BIT)
 357                           ze = - 2 * MANT_BIT - 1;
 358                         else if (ze == 2 * MANT_BIT)
 359                           ze = MAX_EXP - MANT_BIT - 1;
 360                       }
 361                 }
 362           }
 363   }
 364   /* TODO: Similar tests with
 365      x = (-1)^xs * (2^xe - 2^(xe-i)), y = (-1)^ys * (2^ye - 2^(ye-i))  */
 366   /* A product x * y that consists of one segment of bits (or, if you prefer,
 367      two bits, one with positive weight and one with negative weight).  */
 368   {
 369     volatile DOUBLE x;
 370     volatile DOUBLE y;
 371     volatile DOUBLE z;
 372     volatile DOUBLE result;
 373     volatile DOUBLE expected;
 374     int i;
 375     int xs;
 376     int xe;
 377     int ys;
 378     int ye;
 379     int ze;
 380     DOUBLE sign;
 381
 382     for (i = 1; i <= MANT_BIT - 1; i++)
 383       for (xs = 0; xs < 2; xs++)
 384         for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
 385           {
 386             x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
 387               pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
 388
 389             for (ys = 0; ys < 2; ys++)
 390               for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
 391                 {
 392                   y = /* (-1)^ys * (2^ye - 2^(ye-i)) */
 393                     pow_m1[ys] * (POW2 (ye) - POW2 (ye - i));
 394
 395                   sign = pow_m1[xs + ys];
 396
 397                   /* The exact value of x * y is
 398                      (-1)^(xs+ys) * (2^(xe+ye) - 2^(xe+ye-2*i)) */
 399
 400                   /* Test addition (same signs).  */
 401                   for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
 402                     {
 403                       z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
 404                       result = my_fma (x, y, z);
 405                       if (FORGIVE_DOUBLEDOUBLE_BUG)
 406                         if ((xe + ye == ze + MANT_BIT && i > DBL_MANT_BIT)
 407                             || (xe + ye < ze + MANT_BIT
 408                                 && xe + ye >= ze
 409                                 && i == DBL_MANT_BIT)
 410                             || (xe + ye < ze
 411                                 && xe + ye == ze - MANT_BIT + 2 * i))
 412                           goto skip3;
 413                       if (xe + ye > ze + MANT_BIT + 1)
 414                         {
 415                           if (2 * i > MANT_BIT)
 416                             expected = sign * POW2 (xe + ye);
 417                           else
 418                             expected =
 419                               sign * (POW2 (xe + ye)
 420                                       - POW2 (xe + ye - 2 * i));
 421                         }
 422                       else if (xe + ye == ze + MANT_BIT + 1)
 423                         {
 424                           if (2 * i >= MANT_BIT)
 425                             expected = sign * POW2 (xe + ye);
 426                           else
 427                             expected =
 428                               sign * (POW2 (xe + ye)
 429                                       - POW2 (xe + ye - 2 * i));
 430                         }
 431                       else if (xe + ye >= ze - MANT_BIT + 2 * i)
 432                         {
 433                           if (2 * i > MANT_BIT)
 434                             expected =
 435                               sign * (POW2 (xe + ye) + POW2 (ze));
 436                           else
 437                             expected =
 438                               sign * (POW2 (xe + ye)
 439                                       - POW2 (xe + ye - 2 * i)
 440                                       + POW2 (ze));
 441                         }
 442                       else if (xe + ye >= ze - MANT_BIT + 1)
 443                         expected =
 444                           sign * (POW2 (ze) + POW2 (xe + ye));
 445                       else
 446                         expected = z;
 447                       ASSERT (result == expected);
 448
 449                      skip3:
 450                       ze++;
 451                       /* Shortcut some values of ze, to speed up the test.  */
 452                       if (ze == MIN_EXP + MANT_BIT)
 453                         ze = - 2 * MANT_BIT - 1;
 454                       else if (ze == 2 * MANT_BIT)
 455                         ze = MAX_EXP - MANT_BIT - 1;
 456                     }
 457
 458                   /* Test subtraction (opposite signs).  */
 459                   if (i > 1)
 460                     for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
 461                       {
 462                         z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
 463                         result = my_fma (x, y, z);
 464                         if (FORGIVE_DOUBLEDOUBLE_BUG)
 465                           if (xe + ye > ze
 466                               && xe + ye < ze + DBL_MANT_BIT
 467                               && xe + ye == ze + 2 * i - LDBL_MANT_BIT)
 468                             goto skip4;
 469                         if (xe + ye == ze)
 470                           {
 471                             /* maximal extinction */
 472                             expected = sign * - POW2 (xe + ye - 2 * i);
 473                           }
 474                         else if (xe + ye > ze + MANT_BIT + 1)
 475                           {
 476                             if (2 * i > MANT_BIT + 1)
 477                               expected = sign * POW2 (xe + ye);
 478                             else if (2 * i == MANT_BIT + 1)
 479                               expected =
 480                                 sign * (POW2 (xe + ye)
 481                                         - POW2 (xe + ye - MANT_BIT));
 482                             else
 483                               expected =
 484                                 sign * (POW2 (xe + ye)
 485                                         - POW2 (xe + ye - 2 * i));
 486                           }
 487                         else if (xe + ye == ze + MANT_BIT + 1)
 488                           {
 489                             if (2 * i > MANT_BIT)
 490                               expected =
 491                                 sign * (POW2 (xe + ye)
 492                                         - POW2 (xe + ye - MANT_BIT));
 493                             else if (2 * i == MANT_BIT)
 494                               expected =
 495                                 sign * (POW2 (xe + ye)
 496                                         - POW2 (xe + ye - MANT_BIT + 1));
 497                             else
 498                               expected =
 499                                 sign * (POW2 (xe + ye)
 500                                         - POW2 (xe + ye - 2 * i));
 501                           }
 502                         else if (xe + ye == ze + MANT_BIT)
 503                           {
 504                             if (2 * i > MANT_BIT + 1)
 505                               expected =
 506                                 sign * (POW2 (xe + ye)
 507                                         - POW2 (xe + ye - MANT_BIT));
 508                             else if (2 * i == MANT_BIT + 1)
 509                               expected =
 510                                 sign * (POW2 (xe + ye)
 511                                         - POW2 (xe + ye - MANT_BIT + 1));
 512                             else
 513                               expected =
 514                                 sign * (POW2 (xe + ye)
 515                                         - POW2 (ze)
 516                                         - POW2 (xe + ye - 2 * i));
 517                           }
 518                         else if (xe + ye > ze - MANT_BIT + 2 * i)
 519                           {
 520                             if (2 * i > MANT_BIT)
 521                               expected =
 522                                 sign * (POW2 (xe + ye) - POW2 (ze));
 523                             else
 524                               expected =
 525                                 sign * (POW2 (xe + ye)
 526                                         - POW2 (ze)
 527                                         - POW2 (xe + ye - 2 * i));
 528                           }
 529                         else if (xe + ye == ze - MANT_BIT + 2 * i)
 530                           expected =
 531                             sign * (- POW2 (ze)
 532                                     + POW2 (xe + ye)
 533                                     - POW2 (xe + ye - 2 * i));
 534                         else if (xe + ye >= ze - MANT_BIT)
 535                           expected = sign * (- POW2 (ze) + POW2 (xe + ye));
 536                         else
 537                           expected = z;
 538                         ASSERT (result == expected);
 539
 540                        skip4:
 541                         ze++;
 542                         /* Shortcut some values of ze, to speed up the test.  */
 543                         if (ze == MIN_EXP + MANT_BIT)
 544                           ze = - 2 * MANT_BIT - 1;
 545                         else if (ze == 2 * MANT_BIT)
 546                           ze = MAX_EXP - MANT_BIT - 1;
 547                       }
 548                 }
 549           }
 550   }
 551   /* TODO: Tests with denormalized results.  */
 552   /* Tests with temporary overflow.  */
 553   {
 554     volatile DOUBLE x = POW2 (MAX_EXP - 1);
 555     volatile DOUBLE y = POW2 (MAX_EXP - 1);
 556     volatile DOUBLE z = - INFINITY;
 557     volatile DOUBLE result = my_fma (x, y, z);
 558     ASSERT (result == - INFINITY);
 559   }
 560   {
 561     volatile DOUBLE x = POW2 (MAX_EXP - 1); /* 2^(MAX_EXP-1) */
 562     volatile DOUBLE y = L_(2.0);            /* 2^1 */
 563     volatile DOUBLE z =               /* -(2^MAX_EXP - 2^(MAX_EXP-MANT_BIT)) */
 564       - LDEXP (POW2 (MAX_EXP - 1) - POW2 (MAX_EXP - MANT_BIT - 1), 1);
 565     volatile DOUBLE result = my_fma (x, y, z);
 566     if (!FORGIVE_DOUBLEDOUBLE_BUG)
 567       ASSERT (result == POW2 (MAX_EXP - MANT_BIT));
 568   }
 569   {
 570     volatile DOUBLE x = POW2 (MAX_EXP - 1);             /* 2^(MAX_EXP-1) */
 571     volatile DOUBLE y = L_(3.0);                        /* 3 */
 572     volatile DOUBLE z = - LDEXP (L_(5.0), MAX_EXP - 3); /* -5*2^(MAX_EXP-3) */
 573     volatile DOUBLE result = my_fma (x, y, z);
 574     ASSERT (result == LDEXP (L_(7.0), MAX_EXP - 3));
 575   }
 576 }