This source file includes following definitions.
- test_function
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 #define POW2(e) \
21 LDEXP (L_(1.0), e)
22
23
24
25 #define XE_RANGE 0
26 #define YE_RANGE 0
27
28
29
30
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
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
45
46 static void
47 test_function (DOUBLE (*my_fma) (DOUBLE, DOUBLE, DOUBLE))
48 {
49
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
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);
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);
79
80 sign = pow_m1[xs + ys];
81
82
83 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
84 {
85 z = sign * POW2 (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
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
104 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
105 {
106 z = - sign * POW2 (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
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
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 =
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 =
154 pow_m1[ys] * (POW2 (ye) + POW2 (ye - i));
155
156 sign = pow_m1[xs + ys];
157
158
159
160
161
162 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
163 {
164 z = sign * POW2 (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
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
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
258 if (i > 1)
259 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
260 {
261 z = - sign * POW2 (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
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
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
315 expected =
316 sign * (POW2 (xe + ye)
317 + POW2 (xe + ye - i + 1));
318 else
319
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
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
365
366
367
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 =
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 =
393 pow_m1[ys] * (POW2 (ye) - POW2 (ye - i));
394
395 sign = pow_m1[xs + ys];
396
397
398
399
400
401 for (ze = MIN_EXP - MANT_BIT; ze <= MAX_EXP - 1;)
402 {
403 z = sign * POW2 (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
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
459 if (i > 1)
460 for (ze = MIN_SUBNORMAL_EXP; ze <= MAX_EXP - 1;)
461 {
462 z = - sign * POW2 (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
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
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
552
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);
562 volatile DOUBLE y = L_(2.0);
563 volatile DOUBLE z =
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);
571 volatile DOUBLE y = L_(3.0);
572 volatile DOUBLE z = - LDEXP (L_(5.0), MAX_EXP - 3);
573 volatile DOUBLE result = my_fma (x, y, z);
574 ASSERT (result == LDEXP (L_(7.0), MAX_EXP - 3));
575 }
576 }