1/* Test floating-point arithmetic operations.
2 Copyright (C) 1997-2022 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
9
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <https://www.gnu.org/licenses/>. */
18#ifndef _GNU_SOURCE
19#define _GNU_SOURCE
20#endif
21#include <math.h>
22#include <stdio.h>
23#include <stdlib.h>
24#include <string.h>
25#include <fenv.h>
26#include <assert.h>
27
28#ifndef ESIZE
29typedef double tocheck_t;
30#define ESIZE 11
31#define MSIZE 52
32#define FUNC(x) x
33#endif
34
35#define R_NEAREST 1
36#define R_ZERO 2
37#define R_UP 4
38#define R_DOWN 8
39#define R_ALL (R_NEAREST|R_ZERO|R_UP|R_DOWN)
40static fenv_t rmodes[4];
41static const char * const rmnames[4] =
42{ "nearest","zero","+Inf","-Inf" };
43
44typedef union {
45 tocheck_t tc;
46 unsigned char c[sizeof (tocheck_t)];
47} union_t;
48
49/* Don't try reading these in a font that doesn't distinguish
50 O and zero. */
51typedef enum {
52 P_Z = 0x0, /* 00000...0 */
53 P_000O = 0x1, /* 00011...1 */
54 P_001Z = 0x2, /* 00100...0 */
55 P_00O = 0x3, /* 00111...1 */
56 P_01Z = 0x4, /* 01000...0 */
57 P_010O = 0x5, /* 01011...1 */
58 P_011Z = 0x6, /* 01100...0 */
59 P_0O = 0x7, /* 01111...1 */
60 P_1Z = 0x8, /* 10000...0 */
61 P_100O = 0x9, /* 10011...1 */
62 P_101Z = 0xa, /* 10100...0 */
63 P_10O = 0xb, /* 10111...1 */
64 P_11Z = 0xc, /* 11000...0 */
65 P_110O = 0xd, /* 11011...1 */
66 P_111Z = 0xe, /* 11100...0 */
67 P_O = 0xf, /* 11111...1 */
68 P_Z1 = 0x11, /* 000...001 */
69 P_Z10 = 0x12, /* 000...010 */
70 P_Z11 = 0x13, /* 000...011 */
71 P_0O00 = 0x14, /* 011...100 */
72 P_0O01 = 0x15, /* 011...101 */
73 P_0O0 = 0x16, /* 011...110 */
74 P_1Z1 = 0x19, /* 100...001 */
75 P_1Z10 = 0x1a, /* 100...010 */
76 P_1Z11 = 0x1b, /* 100...011 */
77 P_O00 = 0x1c, /* 111...100 */
78 P_O01 = 0x1d, /* 111...101 */
79 P_O0 = 0x1e, /* 111...110 */
80 P_R = 0x20, /* rrr...rrr */ /* ('r' means random. ) */
81 P_Ro = 0x21, /* rrr...rrr, with odd parity. */
82 P_0R = 0x22, /* 0rr...rrr */
83 P_1R = 0x23, /* 1rr...rrr */
84 P_Rno = 0x24, /* rrr...rrr, but not all ones. */
85} pattern_t;
86
87static void
88pattern_fill(pattern_t ptn, unsigned char *start, int bitoffset, int count)
89{
90#define bitset(count, value) \
91 start[(count)/8] = (start[(count)/8] & ~(1 << 7-(count)%8) \
92 | (value) << 7-(count)%8)
93 int i;
94
95 if (ptn >= 0 && ptn <= 0xf)
96 {
97 /* Patterns between 0 and 0xF have the following format:
98 The LSBit is used to fill the last n-3 bits of the pattern;
99 The next 3 bits are the first 3 bits of the pattern. */
100 for (i = 0; i < count; i++)
101 if (i < 3)
102 bitset((bitoffset+i), ptn >> (3-i) & 1);
103 else
104 bitset((bitoffset+i), ptn >> 0 & 1);
105 }
106 else if (ptn <= 0x1f)
107 {
108 /* Patterns between 0x10 and 0x1F have the following format:
109 The two LSBits are the last two bits of the pattern;
110 The 0x8 bit is the first bit of the pattern;
111 The 0x4 bit is used to fill the remainder. */
112 for (i = 0; i < count; i++)
113 if (i == 0)
114 bitset((bitoffset+i), ptn >> 3 & 1);
115 else if (i >= count-2)
116 bitset((bitoffset+i), ptn >> (count-1-i) & 1);
117 else
118 bitset((bitoffset+i), ptn >> 2 & 1);
119 }
120 else switch (ptn)
121 {
122 case P_0R: case P_1R:
123 assert(count > 0);
124 bitset(bitoffset, ptn & 1);
125 count--;
126 bitoffset++;
127 case P_R:
128 for (; count > 0; count--, bitoffset++)
129 bitset(bitoffset, rand() & 1);
130 break;
131 case P_Ro:
132 {
133 int op = 1;
134 assert(count > 0);
135 for (; count > 1; count--, bitoffset++)
136 bitset(bitoffset, op ^= (rand() & 1));
137 bitset(bitoffset, op);
138 break;
139 }
140 case P_Rno:
141 {
142 int op = 1;
143 assert(count > 0);
144 for (; count > 1; count--, bitoffset++)
145 {
146 int r = rand() & 1;
147 op &= r;
148 bitset(bitoffset, r);
149 }
150 bitset(bitoffset, rand() & (op ^ 1));
151 break;
152 }
153
154 default:
155 assert(0);
156 }
157#undef bitset
158}
159
160static tocheck_t
161pattern(int negative, pattern_t exp, pattern_t mant)
162{
163 union_t result;
164#if 0
165 int i;
166#endif
167
168 pattern_fill(ptn: negative ? P_O : P_Z, start: result.c, bitoffset: 0, count: 1);
169 pattern_fill(ptn: exp, start: result.c, bitoffset: 1, ESIZE);
170 pattern_fill(ptn: mant, start: result.c, ESIZE+1, MSIZE);
171#if 0
172 printf("neg=%d exp=%02x mant=%02x: ", negative, exp, mant);
173 for (i = 0; i < sizeof (tocheck_t); i++)
174 printf("%02x", result.c[i]);
175 printf("\n");
176#endif
177 return result.tc;
178}
179
180/* Return the closest different tocheck_t to 'x' in the direction of
181 'direction', or 'x' if there is no such value. Assumes 'x' is not
182 a NaN. */
183static tocheck_t
184delta(tocheck_t x, int direction)
185{
186 union_t xx;
187 int i;
188
189 xx.tc = x;
190 if (xx.c[0] & 0x80)
191 direction = -direction;
192 if (direction == +1)
193 {
194 union_t tx;
195 tx.tc = pattern(negative: xx.c[0] >> 7, exp: P_O, mant: P_Z);
196 if (memcmp (tx.c, xx.c, sizeof (tocheck_t)) == 0)
197 return x;
198 }
199 for (i = sizeof (tocheck_t)-1; i > 0; i--)
200 {
201 xx.c[i] += direction;
202 if (xx.c[i] != (direction > 0 ? 0 : 0xff))
203 return xx.tc;
204 }
205 if (direction < 0 && (xx.c[0] & 0x7f) == 0)
206 return pattern(negative: ~(xx.c[0] >> 7) & 1, exp: P_Z, mant: P_Z1);
207 else
208 {
209 xx.c[0] += direction;
210 return xx.tc;
211 }
212}
213
214static int nerrors = 0;
215
216#ifdef FE_ALL_INVALID
217static const int all_exceptions = FE_ALL_INVALID | FE_ALL_EXCEPT;
218#else
219static const int all_exceptions = FE_ALL_EXCEPT;
220#endif
221
222static void
223check_result(int line, const char *rm, tocheck_t expected, tocheck_t actual)
224{
225 if (memcmp (&expected, &actual, sizeof (tocheck_t)) != 0)
226 {
227 unsigned char *ex, *ac;
228 size_t i;
229
230 printf(format: "%s:%d:round %s:result failed\n"
231 " expected result 0x", __FILE__, line, rm);
232 ex = (unsigned char *)&expected;
233 ac = (unsigned char *)&actual;
234 for (i = 0; i < sizeof (tocheck_t); i++)
235 printf(format: "%02x", ex[i]);
236 printf(format: " got 0x");
237 for (i = 0; i < sizeof (tocheck_t); i++)
238 printf(format: "%02x", ac[i]);
239 printf(format: "\n");
240 nerrors++;
241 }
242}
243
244static const struct {
245 int except;
246 const char *name;
247} excepts[] = {
248#define except_entry(ex) { ex, #ex } ,
249#ifdef FE_INEXACT
250 except_entry(FE_INEXACT)
251#else
252# define FE_INEXACT 0
253#endif
254#ifdef FE_DIVBYZERO
255 except_entry(FE_DIVBYZERO)
256#else
257# define FE_DIVBYZERO 0
258#endif
259#ifdef FE_UNDERFLOW
260 except_entry(FE_UNDERFLOW)
261#else
262# define FE_UNDERFLOW 0
263#endif
264#ifdef FE_OVERFLOW
265 except_entry(FE_OVERFLOW)
266#else
267# define FE_OVERFLOW 0
268#endif
269#ifdef FE_INVALID
270 except_entry(FE_INVALID)
271#else
272# define FE_INVALID 0
273#endif
274#ifdef FE_INVALID_SNAN
275 except_entry(FE_INVALID_SNAN)
276#else
277# define FE_INVALID_SNAN FE_INVALID
278#endif
279#ifdef FE_INVALID_ISI
280 except_entry(FE_INVALID_ISI)
281#else
282# define FE_INVALID_ISI FE_INVALID
283#endif
284#ifdef FE_INVALID_IDI
285 except_entry(FE_INVALID_IDI)
286#else
287# define FE_INVALID_IDI FE_INVALID
288#endif
289#ifdef FE_INVALID_ZDZ
290 except_entry(FE_INVALID_ZDZ)
291#else
292# define FE_INVALID_ZDZ FE_INVALID
293#endif
294#ifdef FE_INVALID_COMPARE
295 except_entry(FE_INVALID_COMPARE)
296#else
297# define FE_INVALID_COMPARE FE_INVALID
298#endif
299#ifdef FE_INVALID_SOFTWARE
300 except_entry(FE_INVALID_SOFTWARE)
301#else
302# define FE_INVALID_SOFTWARE FE_INVALID
303#endif
304#ifdef FE_INVALID_SQRT
305 except_entry(FE_INVALID_SQRT)
306#else
307# define FE_INVALID_SQRT FE_INVALID
308#endif
309#ifdef FE_INVALID_INTEGER_CONVERSION
310 except_entry(FE_INVALID_INTEGER_CONVERSION)
311#else
312# define FE_INVALID_INTEGER_CONVERSION FE_INVALID
313#endif
314};
315
316static int excepts_missing = 0;
317
318static void
319check_excepts(int line, const char *rm, int expected, int actual)
320{
321 if (expected & excepts_missing)
322 expected = expected & ~excepts_missing | FE_INVALID_SNAN;
323 if ((expected & all_exceptions) != actual)
324 {
325 size_t i;
326 printf(format: "%s:%d:round %s:exceptions failed\n"
327 " expected exceptions ", __FILE__, line,rm);
328 for (i = 0; i < sizeof (excepts) / sizeof (excepts[0]); i++)
329 if (expected & excepts[i].except)
330 printf(format: "%s ",excepts[i].name);
331 if ((expected & all_exceptions) == 0)
332 printf(format: "- ");
333 printf(format: "got");
334 for (i = 0; i < sizeof (excepts) / sizeof (excepts[0]); i++)
335 if (actual & excepts[i].except)
336 printf(format: " %s",excepts[i].name);
337 if ((actual & all_exceptions) == 0)
338 printf(format: "- ");
339 printf(format: ".\n");
340 nerrors++;
341 }
342}
343
344typedef enum {
345 B_ADD, B_SUB, B_MUL, B_DIV, B_NEG, B_ABS, B_SQRT
346} op_t;
347typedef struct {
348 int line;
349 op_t op;
350 int a_sgn;
351 pattern_t a_exp, a_mant;
352 int b_sgn;
353 pattern_t b_exp, b_mant;
354 int rmode;
355 int excepts;
356 int x_sgn;
357 pattern_t x_exp, x_mant;
358} optest_t;
359static const optest_t optests[] = {
360 /* Additions of zero. */
361 {__LINE__,B_ADD, 0,P_Z,P_Z, 0,P_Z,P_Z, R_ALL,0, 0,P_Z,P_Z },
362 {__LINE__,B_ADD, 1,P_Z,P_Z, 0,P_Z,P_Z, R_ALL & ~R_DOWN,0, 0,P_Z,P_Z },
363 {__LINE__,B_ADD, 1,P_Z,P_Z, 0,P_Z,P_Z, R_DOWN,0, 1,P_Z,P_Z },
364 {__LINE__,B_ADD, 1,P_Z,P_Z, 1,P_Z,P_Z, R_ALL,0, 1,P_Z,P_Z },
365
366 /* Additions with NaN. */
367 {__LINE__,B_ADD, 0,P_O,P_101Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_101Z },
368 {__LINE__,B_ADD, 0,P_O,P_01Z, 0,P_Z,P_Z, R_ALL,
369 FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_11Z },
370 {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_O,P_0O, R_ALL,
371 FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_O },
372 {__LINE__,B_ADD, 0,P_Z,P_Z, 0,P_O,P_11Z, R_ALL,0, 0,P_O,P_11Z },
373 {__LINE__,B_ADD, 0,P_O,P_001Z, 0,P_O,P_001Z, R_ALL,
374 FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_101Z },
375 {__LINE__,B_ADD, 0,P_O,P_1Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_1Z },
376 {__LINE__,B_ADD, 0,P_0O,P_Z, 0,P_O,P_10O, R_ALL,0, 0,P_O,P_10O },
377
378 /* Additions with infinity. */
379 {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_Z },
380 {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_Z,P_Z, R_ALL,0, 0,P_O,P_Z },
381 {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_Z,P_Z, R_ALL,0, 1,P_O,P_Z },
382 {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_Z,P_Z, R_ALL,0, 1,P_O,P_Z },
383 {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_O,P_Z, R_ALL,0, 0,P_O,P_Z },
384 {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_O,P_Z, R_ALL,0, 1,P_O,P_Z },
385 {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_O,P_Z, R_ALL,
386 FE_INVALID | FE_INVALID_ISI, 0,P_O,P_1Z },
387 {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_O,P_Z, R_ALL,
388 FE_INVALID | FE_INVALID_ISI, 0,P_O,P_1Z },
389 {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_0O,P_Z, R_ALL,0, 0,P_O,P_Z },
390 {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_0O,P_Z, R_ALL,0, 1,P_O,P_Z },
391 {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_0O,P_Z, R_ALL,0, 0,P_O,P_Z },
392 {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_0O,P_Z, R_ALL,0, 1,P_O,P_Z },
393
394 /* Overflow (and zero). */
395 {__LINE__,B_ADD, 0,P_O0,P_Z, 0,P_O0,P_Z, R_NEAREST | R_UP,
396 FE_INEXACT | FE_OVERFLOW, 0,P_O,P_Z },
397 {__LINE__,B_ADD, 0,P_O0,P_Z, 0,P_O0,P_Z, R_ZERO | R_DOWN,
398 FE_INEXACT | FE_OVERFLOW, 0,P_O0,P_O },
399 {__LINE__,B_ADD, 1,P_O0,P_Z, 1,P_O0,P_Z, R_NEAREST | R_DOWN,
400 FE_INEXACT | FE_OVERFLOW, 1,P_O,P_Z },
401 {__LINE__,B_ADD, 1,P_O0,P_Z, 1,P_O0,P_Z, R_ZERO | R_UP,
402 FE_INEXACT | FE_OVERFLOW, 1,P_O0,P_O },
403 {__LINE__,B_ADD, 0,P_O0,P_Z, 1,P_O0,P_Z, R_ALL & ~R_DOWN,
404 0, 0,P_Z,P_Z },
405 {__LINE__,B_ADD, 0,P_O0,P_Z, 1,P_O0,P_Z, R_DOWN,
406 0, 1,P_Z,P_Z },
407
408 /* Negation. */
409 {__LINE__,B_NEG, 0,P_Z,P_Z, 0,0,0, R_ALL, 0, 1,P_Z,P_Z },
410 {__LINE__,B_NEG, 1,P_Z,P_Z, 0,0,0, R_ALL, 0, 0,P_Z,P_Z },
411 {__LINE__,B_NEG, 0,P_O,P_Z, 0,0,0, R_ALL, 0, 1,P_O,P_Z },
412 {__LINE__,B_NEG, 1,P_O,P_Z, 0,0,0, R_ALL, 0, 0,P_O,P_Z },
413 {__LINE__,B_NEG, 0,P_O,P_1Z, 0,0,0, R_ALL, 0, 1,P_O,P_1Z },
414 {__LINE__,B_NEG, 1,P_O,P_1Z, 0,0,0, R_ALL, 0, 0,P_O,P_1Z },
415 {__LINE__,B_NEG, 0,P_O,P_01Z, 0,0,0, R_ALL, 0, 1,P_O,P_01Z },
416 {__LINE__,B_NEG, 1,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
417 {__LINE__,B_NEG, 0,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 1,P_1Z,P_1Z1 },
418 {__LINE__,B_NEG, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
419 {__LINE__,B_NEG, 0,P_Z,P_Z1, 0,0,0, R_ALL, 0, 1,P_Z,P_Z1 },
420 {__LINE__,B_NEG, 1,P_Z,P_Z1, 0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
421
422 /* Absolute value. */
423 {__LINE__,B_ABS, 0,P_Z,P_Z, 0,0,0, R_ALL, 0, 0,P_Z,P_Z },
424 {__LINE__,B_ABS, 1,P_Z,P_Z, 0,0,0, R_ALL, 0, 0,P_Z,P_Z },
425 {__LINE__,B_ABS, 0,P_O,P_Z, 0,0,0, R_ALL, 0, 0,P_O,P_Z },
426 {__LINE__,B_ABS, 1,P_O,P_Z, 0,0,0, R_ALL, 0, 0,P_O,P_Z },
427 {__LINE__,B_ABS, 0,P_O,P_1Z, 0,0,0, R_ALL, 0, 0,P_O,P_1Z },
428 {__LINE__,B_ABS, 1,P_O,P_1Z, 0,0,0, R_ALL, 0, 0,P_O,P_1Z },
429 {__LINE__,B_ABS, 0,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
430 {__LINE__,B_ABS, 1,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
431 {__LINE__,B_ABS, 0,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
432 {__LINE__,B_ABS, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
433 {__LINE__,B_ABS, 0,P_Z,P_Z1, 0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
434 {__LINE__,B_ABS, 1,P_Z,P_Z1, 0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
435
436 /* Square root. */
437 {__LINE__,B_SQRT, 0,P_Z,P_Z, 0,0,0, R_ALL, 0, 0,P_Z,P_Z },
438 {__LINE__,B_SQRT, 1,P_Z,P_Z, 0,0,0, R_ALL, 0, 1,P_Z,P_Z },
439 {__LINE__,B_SQRT, 0,P_O,P_1Z, 0,0,0, R_ALL, 0, 0,P_O,P_1Z },
440 {__LINE__,B_SQRT, 1,P_O,P_1Z, 0,0,0, R_ALL, 0, 1,P_O,P_1Z },
441 {__LINE__,B_SQRT, 0,P_O,P_01Z, 0,0,0, R_ALL,
442 FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_11Z },
443 {__LINE__,B_SQRT, 1,P_O,P_01Z, 0,0,0, R_ALL,
444 FE_INVALID | FE_INVALID_SNAN, 1,P_O,P_11Z },
445
446 {__LINE__,B_SQRT, 0,P_O,P_Z, 0,0,0, R_ALL, 0, 0,P_O,P_Z },
447 {__LINE__,B_SQRT, 0,P_0O,P_Z, 0,0,0, R_ALL, 0, 0,P_0O,P_Z },
448
449 {__LINE__,B_SQRT, 1,P_O,P_Z, 0,0,0, R_ALL,
450 FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
451 {__LINE__,B_SQRT, 1,P_1Z,P_1Z1, 0,0,0, R_ALL,
452 FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
453 {__LINE__,B_SQRT, 1,P_Z,P_Z1, 0,0,0, R_ALL,
454 FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
455
456};
457
458static void
459check_op(void)
460{
461 size_t i;
462 int j;
463 tocheck_t r, a, b, x;
464 int raised;
465
466 for (i = 0; i < sizeof (optests) / sizeof (optests[0]); i++)
467 {
468 a = pattern(negative: optests[i].a_sgn, exp: optests[i].a_exp,
469 mant: optests[i].a_mant);
470 b = pattern(negative: optests[i].b_sgn, exp: optests[i].b_exp,
471 mant: optests[i].b_mant);
472 x = pattern(negative: optests[i].x_sgn, exp: optests[i].x_exp,
473 mant: optests[i].x_mant);
474 for (j = 0; j < 4; j++)
475 if (optests[i].rmode & 1<<j)
476 {
477 fesetenv(envp: rmodes+j);
478 switch (optests[i].op)
479 {
480 case B_ADD: r = a + b; break;
481 case B_SUB: r = a - b; break;
482 case B_MUL: r = a * b; break;
483 case B_DIV: r = a / b; break;
484 case B_NEG: r = -a; break;
485 case B_ABS: r = FUNC(fabs)(x: a); break;
486 case B_SQRT: r = FUNC(sqrt)(a); break;
487 }
488 raised = fetestexcept(excepts: all_exceptions);
489 check_result(line: optests[i].line,rm: rmnames[j],expected: x,actual: r);
490 check_excepts(line: optests[i].line,rm: rmnames[j],
491 expected: optests[i].excepts,actual: raised);
492 }
493 }
494}
495
496static void
497fail_xr(int line, const char *rm, tocheck_t x, tocheck_t r, tocheck_t xx,
498 int xflag)
499{
500 size_t i;
501 unsigned char *cx, *cr, *cxx;
502
503 printf(format: "%s:%d:round %s:fail\n with x=0x", __FILE__, line,rm);
504 cx = (unsigned char *)&x;
505 cr = (unsigned char *)&r;
506 cxx = (unsigned char *)&xx;
507 for (i = 0; i < sizeof (tocheck_t); i++)
508 printf(format: "%02x", cx[i]);
509 printf(format: " r=0x");
510 for (i = 0; i < sizeof (tocheck_t); i++)
511 printf(format: "%02x", cr[i]);
512 printf(format: " xx=0x");
513 for (i = 0; i < sizeof (tocheck_t); i++)
514 printf(format: "%02x", cxx[i]);
515 printf(format: " inexact=%d\n", xflag != 0);
516 nerrors++;
517}
518
519static void
520check_sqrt(tocheck_t a)
521{
522 int j;
523 tocheck_t r0, r1, r2, x0, x1, x2;
524 int raised = 0;
525 int ok;
526
527 for (j = 0; j < 4; j++)
528 {
529 int excepts;
530
531 fesetenv(envp: rmodes+j);
532 r1 = FUNC(sqrt)(a);
533 excepts = fetestexcept(excepts: all_exceptions);
534 fesetenv(FE_DFL_ENV);
535 raised |= excepts & ~FE_INEXACT;
536 x1 = r1 * r1 - a;
537 if (excepts & FE_INEXACT)
538 {
539 r0 = delta(x: r1,direction: -1); r2 = delta(x: r1,direction: 1);
540 switch (1 << j)
541 {
542 case R_NEAREST:
543 x0 = r0 * r0 - a; x2 = r2 * r2 - a;
544 ok = fabs(x: x0) >= fabs(x: x1) && fabs(x: x1) <= fabs(x: x2);
545 break;
546 case R_ZERO: case R_DOWN:
547 x2 = r2 * r2 - a;
548 ok = x1 <= 0 && x2 >= 0;
549 break;
550 case R_UP:
551 x0 = r0 * r0 - a;
552 ok = x1 >= 0 && x0 <= 0;
553 break;
554 default:
555 assert(0);
556 }
557 }
558 else
559 ok = x1 == 0;
560 if (!ok)
561 fail_xr(__LINE__,rm: rmnames[j],x: a,r: r1,xx: x1,xflag: excepts&FE_INEXACT);
562 }
563 check_excepts(__LINE__,rm: "all",expected: 0,actual: raised);
564}
565
566int main(int argc, char **argv)
567{
568 int i;
569
570 _LIB_VERSION = _IEEE_;
571
572 /* Set up environments for rounding modes. */
573 fesetenv(FE_DFL_ENV);
574 fesetround(FE_TONEAREST);
575 fegetenv(envp: rmodes+0);
576 fesetround(FE_TOWARDZERO);
577 fegetenv(envp: rmodes+1);
578 fesetround(FE_UPWARD);
579 fegetenv(envp: rmodes+2);
580 fesetround(FE_DOWNWARD);
581 fegetenv(envp: rmodes+3);
582
583#if defined(FE_INVALID_SOFTWARE) || defined(FE_INVALID_SQRT)
584 /* There's this really stupid feature of the 601... */
585 fesetenv(FE_DFL_ENV);
586 feraiseexcept(FE_INVALID_SOFTWARE);
587 if (!fetestexcept(FE_INVALID_SOFTWARE))
588 excepts_missing |= FE_INVALID_SOFTWARE;
589 fesetenv(FE_DFL_ENV);
590 feraiseexcept(FE_INVALID_SQRT);
591 if (!fetestexcept(FE_INVALID_SQRT))
592 excepts_missing |= FE_INVALID_SQRT;
593#endif
594
595 check_op();
596 for (i = 0; i < 100000; i++)
597 check_sqrt(a: pattern(negative: 0, exp: P_Rno, mant: P_R));
598 for (i = 0; i < 100; i++)
599 check_sqrt(a: pattern(negative: 0, exp: P_Z, mant: P_R));
600 check_sqrt(a: pattern(negative: 0,exp: P_Z,mant: P_Z1));
601
602 printf(format: "%d errors.\n", nerrors);
603 return nerrors == 0 ? 0 : 1;
604}
605

source code of glibc/sysdeps/powerpc/test-arith.c