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