1 /* Test of fused multiply-add.
2 Copyright (C) 2011 Free Software Foundation, Inc.
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.
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.
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/>. */
17 /* Written by Bruno Haible <bruno@clisp.org>, 2011. */
19 /* Returns 2^e as a DOUBLE. */
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. */
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
34 # define FORGIVE_DOUBLEDOUBLE_BUG 0
37 /* Subnormal numbers appear to not work as expected on IRIX 6.5. */
39 # define MIN_SUBNORMAL_EXP (MIN_EXP - 1)
41 # define MIN_SUBNORMAL_EXP (MIN_EXP - MANT_BIT)
44 /* Check rounding behaviour. */
47 test_function (DOUBLE (*my_fma) (DOUBLE, DOUBLE, DOUBLE))
49 /* Array mapping n to (-1)^n. */
50 static const DOUBLE pow_m1[] =
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)
58 volatile DOUBLE result;
59 volatile DOUBLE expected;
61 /* A product x * y that consists of two bits. */
70 for (xs = 0; xs < 2; xs++)
71 for (xe = -XE_RANGE; xe <= XE_RANGE; xe++)
73 x = pow_m1[xs] * POW2 (xe); /* (-1)^xs * 2^xe */
75 for (ys = 0; ys < 2; ys++)
76 for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
78 y = pow_m1[ys] * POW2 (ye); /* (-1)^ys * 2^ye */
80 sign = pow_m1[xs + ys];
82 /* Test addition (same signs). */
83 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
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));
93 ASSERT (result == expected);
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;
103 /* Test subtraction (opposite signs). */
104 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
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));
116 ASSERT (result == expected);
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;
128 /* A product x * y that consists of three bits. */
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++)
142 x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
143 pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
145 for (ys = 0; ys < 2; ys++)
146 for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
148 y = /* (-1)^ys * (2^ye + 2^(ye-i)) */
149 pow_m1[ys] * (POW2 (ye) + POW2 (ye - i));
151 sign = pow_m1[xs + ys];
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)) */
156 /* Test addition (same signs). */
157 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
159 z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
160 result = my_fma (x, y, z);
161 if (FORGIVE_DOUBLEDOUBLE_BUG)
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))
168 if (xe + ye > ze + MANT_BIT)
170 if (2 * i > MANT_BIT)
172 sign * (POW2 (xe + ye)
173 + POW2 (xe + ye - i + 1));
174 else if (2 * i == MANT_BIT)
176 sign * (POW2 (xe + ye)
177 + POW2 (xe + ye - i + 1)
178 + POW2 (xe + ye - MANT_BIT + 1));
181 sign * (POW2 (xe + ye)
182 + POW2 (xe + ye - i + 1)
183 + POW2 (xe + ye - 2 * i));
185 else if (xe + ye == ze + MANT_BIT)
187 if (2 * i >= MANT_BIT)
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 */
195 sign * (POW2 (xe + ye)
196 + POW2 (xe + ye - i + 1)
197 + POW2 (xe + ye - 2 * i + 1));
200 sign * (POW2 (xe + ye)
201 + POW2 (xe + ye - i + 1)
202 + POW2 (xe + ye - 2 * i));
204 else if (xe + ye > ze - MANT_BIT + 2 * i)
208 + POW2 (xe + ye - i + 1)
209 + POW2 (xe + ye - 2 * i));
210 else if (xe + ye >= ze - MANT_BIT + i)
214 + POW2 (xe + ye - i + 1));
215 else if (xe + ye == ze - MANT_BIT + i - 1)
219 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
224 + POW2 (ze - MANT_BIT + 1));
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)
230 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
231 else if (xe + ye == ze - MANT_BIT - 1)
235 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
241 ASSERT (result == expected);
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;
252 /* Test subtraction (opposite signs). */
254 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
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)
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))
272 /* maximal extinction */
274 sign * (POW2 (xe + ye - i + 1)
275 + POW2 (xe + ye - 2 * i));
277 else if (xe + ye == ze - 1)
279 /* significant extinction */
280 if (2 * i > MANT_BIT)
282 sign * (- POW2 (xe + ye)
283 + POW2 (xe + ye - i + 1));
286 sign * (- POW2 (xe + ye)
287 + POW2 (xe + ye - i + 1)
288 + POW2 (xe + ye - 2 * i));
290 else if (xe + ye > ze + MANT_BIT)
292 if (2 * i >= MANT_BIT)
294 sign * (POW2 (xe + ye)
295 + POW2 (xe + ye - i + 1));
298 sign * (POW2 (xe + ye)
299 + POW2 (xe + ye - i + 1)
300 + POW2 (xe + ye - 2 * i));
302 else if (xe + ye == ze + MANT_BIT)
304 if (2 * i >= MANT_BIT)
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 */
311 sign * (POW2 (xe + ye)
312 + POW2 (xe + ye - i + 1));
314 /* round-to-even rounds up */
316 sign * (POW2 (xe + ye)
317 + POW2 (xe + ye - i + 1)
318 + POW2 (xe + ye - 2 * i));
320 else if (xe + ye >= ze - MANT_BIT + 2 * i)
324 + POW2 (xe + ye - i + 1)
325 + POW2 (xe + ye - 2 * i));
326 else if (xe + ye >= ze - MANT_BIT + i - 1)
330 + POW2 (xe + ye - i + 1));
331 else if (xe + ye == ze - MANT_BIT + i - 2)
335 + POW2 (ze - MANT_BIT));
336 else if (xe + ye >= ze - MANT_BIT)
340 else if (xe + ye == ze - MANT_BIT - 1)
343 + POW2 (ze - MANT_BIT));
346 ASSERT (result == expected);
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;
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). */
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++)
376 x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
377 pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
379 for (ys = 0; ys < 2; ys++)
380 for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
382 y = /* (-1)^ys * (2^ye - 2^(ye-i)) */
383 pow_m1[ys] * (POW2 (ye) - POW2 (ye - i));
385 sign = pow_m1[xs + ys];
387 /* The exact value of x * y is
388 (-1)^(xs+ys) * (2^(xe+ye) - 2^(xe+ye-2*i)) */
390 /* Test addition (same signs). */
391 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
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
399 && i == DBL_MANT_BIT)
401 && xe + ye == ze - MANT_BIT + 2 * i))
403 if (xe + ye > ze + MANT_BIT + 1)
405 if (2 * i > MANT_BIT)
406 expected = sign * POW2 (xe + ye);
409 sign * (POW2 (xe + ye)
410 - POW2 (xe + ye - 2 * i));
412 else if (xe + ye == ze + MANT_BIT + 1)
414 if (2 * i >= MANT_BIT)
415 expected = sign * POW2 (xe + ye);
418 sign * (POW2 (xe + ye)
419 - POW2 (xe + ye - 2 * i));
421 else if (xe + ye >= ze - MANT_BIT + 2 * i)
423 if (2 * i > MANT_BIT)
425 sign * (POW2 (xe + ye) + POW2 (ze));
428 sign * (POW2 (xe + ye)
429 - POW2 (xe + ye - 2 * i)
432 else if (xe + ye >= ze - MANT_BIT + 1)
434 sign * (POW2 (ze) + POW2 (xe + ye));
437 ASSERT (result == expected);
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;
448 /* Test subtraction (opposite signs). */
450 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
452 z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
453 result = my_fma (x, y, z);
454 if (FORGIVE_DOUBLEDOUBLE_BUG)
456 && xe + ye < ze + DBL_MANT_BIT
457 && xe + ye == ze + 2 * i - LDBL_MANT_BIT)
461 /* maximal extinction */
462 expected = sign * - POW2 (xe + ye - 2 * i);
464 else if (xe + ye > ze + MANT_BIT + 1)
466 if (2 * i > MANT_BIT + 1)
467 expected = sign * POW2 (xe + ye);
468 else if (2 * i == MANT_BIT + 1)
470 sign * (POW2 (xe + ye)
471 - POW2 (xe + ye - MANT_BIT));
474 sign * (POW2 (xe + ye)
475 - POW2 (xe + ye - 2 * i));
477 else if (xe + ye == ze + MANT_BIT + 1)
479 if (2 * i > MANT_BIT)
481 sign * (POW2 (xe + ye)
482 - POW2 (xe + ye - MANT_BIT));
483 else if (2 * i == MANT_BIT)
485 sign * (POW2 (xe + ye)
486 - POW2 (xe + ye - MANT_BIT + 1));
489 sign * (POW2 (xe + ye)
490 - POW2 (xe + ye - 2 * i));
492 else if (xe + ye == ze + MANT_BIT)
494 if (2 * i > MANT_BIT + 1)
496 sign * (POW2 (xe + ye)
497 - POW2 (xe + ye - MANT_BIT));
498 else if (2 * i == MANT_BIT + 1)
500 sign * (POW2 (xe + ye)
501 - POW2 (xe + ye - MANT_BIT + 1));
504 sign * (POW2 (xe + ye)
506 - POW2 (xe + ye - 2 * i));
508 else if (xe + ye > ze - MANT_BIT + 2 * i)
510 if (2 * i > MANT_BIT)
512 sign * (POW2 (xe + ye) - POW2 (ze));
515 sign * (POW2 (xe + ye)
517 - POW2 (xe + ye - 2 * i));
519 else if (xe + ye == ze - MANT_BIT + 2 * i)
523 - POW2 (xe + ye - 2 * i));
524 else if (xe + ye >= ze - MANT_BIT)
525 expected = sign * (- POW2 (ze) + POW2 (xe + ye));
528 ASSERT (result == expected);
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;
541 /* TODO: Tests with denormalized results. */
542 /* Tests with temporary overflow. */
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);
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));
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));