1 /* Test of fused multiply-add.
2 Copyright (C) 2011-2013 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)
56 /* A product x * y that consists of two bits. */
61 volatile DOUBLE result;
62 volatile DOUBLE expected;
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. */
133 volatile DOUBLE result;
134 volatile DOUBLE expected;
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++)
147 x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
148 pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
150 for (ys = 0; ys < 2; ys++)
151 for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
153 y = /* (-1)^ys * (2^ye + 2^(ye-i)) */
154 pow_m1[ys] * (POW2 (ye) + POW2 (ye - i));
156 sign = pow_m1[xs + ys];
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)) */
161 /* Test addition (same signs). */
162 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
164 z = sign * POW2 (ze); /* (-1)^(xs+ys) * 2^ze */
165 result = my_fma (x, y, z);
166 if (FORGIVE_DOUBLEDOUBLE_BUG)
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))
173 if (xe + ye > ze + MANT_BIT)
175 if (2 * i > MANT_BIT)
177 sign * (POW2 (xe + ye)
178 + POW2 (xe + ye - i + 1));
179 else if (2 * i == MANT_BIT)
181 sign * (POW2 (xe + ye)
182 + POW2 (xe + ye - i + 1)
183 + POW2 (xe + ye - MANT_BIT + 1));
186 sign * (POW2 (xe + ye)
187 + POW2 (xe + ye - i + 1)
188 + POW2 (xe + ye - 2 * i));
190 else if (xe + ye == ze + MANT_BIT)
192 if (2 * i >= MANT_BIT)
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 */
200 sign * (POW2 (xe + ye)
201 + POW2 (xe + ye - i + 1)
202 + POW2 (xe + ye - 2 * i + 1));
205 sign * (POW2 (xe + ye)
206 + POW2 (xe + ye - i + 1)
207 + POW2 (xe + ye - 2 * i));
209 else if (xe + ye > ze - MANT_BIT + 2 * i)
213 + POW2 (xe + ye - i + 1)
214 + POW2 (xe + ye - 2 * i));
215 else if (xe + ye >= ze - MANT_BIT + i)
219 + POW2 (xe + ye - i + 1));
220 else if (xe + ye == ze - MANT_BIT + i - 1)
224 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
229 + POW2 (ze - MANT_BIT + 1));
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)
235 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
236 else if (xe + ye == ze - MANT_BIT - 1)
240 sign * (POW2 (ze) + POW2 (ze - MANT_BIT + 1));
246 ASSERT (result == expected);
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;
257 /* Test subtraction (opposite signs). */
259 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
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)
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))
277 /* maximal extinction */
279 sign * (POW2 (xe + ye - i + 1)
280 + POW2 (xe + ye - 2 * i));
282 else if (xe + ye == ze - 1)
284 /* significant extinction */
285 if (2 * i > MANT_BIT)
287 sign * (- POW2 (xe + ye)
288 + POW2 (xe + ye - i + 1));
291 sign * (- POW2 (xe + ye)
292 + POW2 (xe + ye - i + 1)
293 + POW2 (xe + ye - 2 * i));
295 else if (xe + ye > ze + MANT_BIT)
297 if (2 * i >= MANT_BIT)
299 sign * (POW2 (xe + ye)
300 + POW2 (xe + ye - i + 1));
303 sign * (POW2 (xe + ye)
304 + POW2 (xe + ye - i + 1)
305 + POW2 (xe + ye - 2 * i));
307 else if (xe + ye == ze + MANT_BIT)
309 if (2 * i >= MANT_BIT)
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 */
316 sign * (POW2 (xe + ye)
317 + POW2 (xe + ye - i + 1));
319 /* round-to-even rounds up */
321 sign * (POW2 (xe + ye)
322 + POW2 (xe + ye - i + 1)
323 + POW2 (xe + ye - 2 * i));
325 else if (xe + ye >= ze - MANT_BIT + 2 * i)
329 + POW2 (xe + ye - i + 1)
330 + POW2 (xe + ye - 2 * i));
331 else if (xe + ye >= ze - MANT_BIT + i - 1)
335 + POW2 (xe + ye - i + 1));
336 else if (xe + ye == ze - MANT_BIT + i - 2)
340 + POW2 (ze - MANT_BIT));
341 else if (xe + ye >= ze - MANT_BIT)
345 else if (xe + ye == ze - MANT_BIT - 1)
348 + POW2 (ze - MANT_BIT));
351 ASSERT (result == expected);
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;
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). */
372 volatile DOUBLE result;
373 volatile DOUBLE expected;
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++)
386 x = /* (-1)^xs * (2^xe + 2^(xe-i)) */
387 pow_m1[xs] * (POW2 (xe) + POW2 (xe - i));
389 for (ys = 0; ys < 2; ys++)
390 for (ye = -YE_RANGE; ye <= YE_RANGE; ye++)
392 y = /* (-1)^ys * (2^ye - 2^(ye-i)) */
393 pow_m1[ys] * (POW2 (ye) - POW2 (ye - i));
395 sign = pow_m1[xs + ys];
397 /* The exact value of x * y is
398 (-1)^(xs+ys) * (2^(xe+ye) - 2^(xe+ye-2*i)) */
400 /* Test addition (same signs). */
401 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
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
409 && i == DBL_MANT_BIT)
411 && xe + ye == ze - MANT_BIT + 2 * i))
413 if (xe + ye > ze + MANT_BIT + 1)
415 if (2 * i > MANT_BIT)
416 expected = sign * POW2 (xe + ye);
419 sign * (POW2 (xe + ye)
420 - POW2 (xe + ye - 2 * i));
422 else if (xe + ye == ze + MANT_BIT + 1)
424 if (2 * i >= MANT_BIT)
425 expected = sign * POW2 (xe + ye);
428 sign * (POW2 (xe + ye)
429 - POW2 (xe + ye - 2 * i));
431 else if (xe + ye >= ze - MANT_BIT + 2 * i)
433 if (2 * i > MANT_BIT)
435 sign * (POW2 (xe + ye) + POW2 (ze));
438 sign * (POW2 (xe + ye)
439 - POW2 (xe + ye - 2 * i)
442 else if (xe + ye >= ze - MANT_BIT + 1)
444 sign * (POW2 (ze) + POW2 (xe + ye));
447 ASSERT (result == expected);
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;
458 /* Test subtraction (opposite signs). */
460 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
462 z = - sign * POW2 (ze); /* (-1)^(xs+ys+1) * 2^ze */
463 result = my_fma (x, y, z);
464 if (FORGIVE_DOUBLEDOUBLE_BUG)
466 && xe + ye < ze + DBL_MANT_BIT
467 && xe + ye == ze + 2 * i - LDBL_MANT_BIT)
471 /* maximal extinction */
472 expected = sign * - POW2 (xe + ye - 2 * i);
474 else if (xe + ye > ze + MANT_BIT + 1)
476 if (2 * i > MANT_BIT + 1)
477 expected = sign * POW2 (xe + ye);
478 else if (2 * i == MANT_BIT + 1)
480 sign * (POW2 (xe + ye)
481 - POW2 (xe + ye - MANT_BIT));
484 sign * (POW2 (xe + ye)
485 - POW2 (xe + ye - 2 * i));
487 else if (xe + ye == ze + MANT_BIT + 1)
489 if (2 * i > MANT_BIT)
491 sign * (POW2 (xe + ye)
492 - POW2 (xe + ye - MANT_BIT));
493 else if (2 * i == MANT_BIT)
495 sign * (POW2 (xe + ye)
496 - POW2 (xe + ye - MANT_BIT + 1));
499 sign * (POW2 (xe + ye)
500 - POW2 (xe + ye - 2 * i));
502 else if (xe + ye == ze + MANT_BIT)
504 if (2 * i > MANT_BIT + 1)
506 sign * (POW2 (xe + ye)
507 - POW2 (xe + ye - MANT_BIT));
508 else if (2 * i == MANT_BIT + 1)
510 sign * (POW2 (xe + ye)
511 - POW2 (xe + ye - MANT_BIT + 1));
514 sign * (POW2 (xe + ye)
516 - POW2 (xe + ye - 2 * i));
518 else if (xe + ye > ze - MANT_BIT + 2 * i)
520 if (2 * i > MANT_BIT)
522 sign * (POW2 (xe + ye) - POW2 (ze));
525 sign * (POW2 (xe + ye)
527 - POW2 (xe + ye - 2 * i));
529 else if (xe + ye == ze - MANT_BIT + 2 * i)
533 - POW2 (xe + ye - 2 * i));
534 else if (xe + ye >= ze - MANT_BIT)
535 expected = sign * (- POW2 (ze) + POW2 (xe + ye));
538 ASSERT (result == expected);
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;
551 /* TODO: Tests with denormalized results. */
552 /* Tests with temporary overflow. */
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);
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));
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));