![[first]](../icons/first.png){kind=icon}
![[previous]](../icons/n_left.png){kind=icon}
![[next]](../icons/n_right.png){kind=icon}
![[last]](../icons/n_last.png){kind=icon}
![[top]](../icons/n_top.png){kind=icon}
![[bottom]](../icons/bottom.png){kind=icon}
![[index]](../icons/index.png){kind=icon}
![[help]](../icons/help.png){kind=icon}
*/
1 /* Test of fused multiply-add.
2 Copyright (C) 2011-2021 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 <https://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))
/* ![[first]](../icons/n_first.png)
![[top]](../icons/top.png)
*/
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 }
![[last]](../icons/last.png){kind=icon}
![[bottom]](../icons/n_bottom.png){kind=icon}
*/