1/*
2 * semi.c: test implementations of mathlib seminumerical functions
3 *
4 * Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
5 * See https://llvm.org/LICENSE.txt for license information.
6 * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
7 */
8
9#include <stdio.h>
10#include "semi.h"
11
12static void test_rint(uint32 *in, uint32 *out,
13 int isfloor, int isceil) {
14 int sign = in[0] & 0x80000000;
15 int roundup = (isfloor && sign) || (isceil && !sign);
16 uint32 xh, xl, roundword;
17 int ex = (in[0] >> 20) & 0x7FF; /* exponent */
18 int i;
19
20 if ((ex > 0x3ff + 52 - 1) || /* things this big can't be fractional */
21 ((in[0] & 0x7FFFFFFF) == 0 && in[1] == 0)) { /* zero */
22 /* NaN, Inf, a large integer, or zero: just return the input */
23 out[0] = in[0];
24 out[1] = in[1];
25 return;
26 }
27
28 /*
29 * Special case: ex < 0x3ff, ie our number is in (0,1). Return
30 * 1 or 0 according to roundup.
31 */
32 if (ex < 0x3ff) {
33 out[0] = sign | (roundup ? 0x3FF00000 : 0);
34 out[1] = 0;
35 return;
36 }
37
38 /*
39 * We're not short of time here, so we'll do this the hideously
40 * inefficient way. Shift bit by bit so that the units place is
41 * somewhere predictable, round, and shift back again.
42 */
43 xh = in[0];
44 xl = in[1];
45 roundword = 0;
46 for (i = ex; i < 0x3ff + 52; i++) {
47 if (roundword & 1)
48 roundword |= 2; /* preserve sticky bit */
49 roundword = (roundword >> 1) | ((xl & 1) << 31);
50 xl = (xl >> 1) | ((xh & 1) << 31);
51 xh = xh >> 1;
52 }
53 if (roundword && roundup) {
54 xl++;
55 xh += (xl==0);
56 }
57 for (i = ex; i < 0x3ff + 52; i++) {
58 xh = (xh << 1) | ((xl >> 31) & 1);
59 xl = (xl & 0x7FFFFFFF) << 1;
60 }
61 out[0] = xh;
62 out[1] = xl;
63}
64
65char *test_ceil(uint32 *in, uint32 *out) {
66 test_rint(in, out, isfloor: 0, isceil: 1);
67 return NULL;
68}
69
70char *test_floor(uint32 *in, uint32 *out) {
71 test_rint(in, out, isfloor: 1, isceil: 0);
72 return NULL;
73}
74
75static void test_rintf(uint32 *in, uint32 *out,
76 int isfloor, int isceil) {
77 int sign = *in & 0x80000000;
78 int roundup = (isfloor && sign) || (isceil && !sign);
79 uint32 x, roundword;
80 int ex = (*in >> 23) & 0xFF; /* exponent */
81 int i;
82
83 if ((ex > 0x7f + 23 - 1) || /* things this big can't be fractional */
84 (*in & 0x7FFFFFFF) == 0) { /* zero */
85 /* NaN, Inf, a large integer, or zero: just return the input */
86 *out = *in;
87 return;
88 }
89
90 /*
91 * Special case: ex < 0x7f, ie our number is in (0,1). Return
92 * 1 or 0 according to roundup.
93 */
94 if (ex < 0x7f) {
95 *out = sign | (roundup ? 0x3F800000 : 0);
96 return;
97 }
98
99 /*
100 * We're not short of time here, so we'll do this the hideously
101 * inefficient way. Shift bit by bit so that the units place is
102 * somewhere predictable, round, and shift back again.
103 */
104 x = *in;
105 roundword = 0;
106 for (i = ex; i < 0x7F + 23; i++) {
107 if (roundword & 1)
108 roundword |= 2; /* preserve sticky bit */
109 roundword = (roundword >> 1) | ((x & 1) << 31);
110 x = x >> 1;
111 }
112 if (roundword && roundup) {
113 x++;
114 }
115 for (i = ex; i < 0x7F + 23; i++) {
116 x = x << 1;
117 }
118 *out = x;
119}
120
121char *test_ceilf(uint32 *in, uint32 *out) {
122 test_rintf(in, out, isfloor: 0, isceil: 1);
123 return NULL;
124}
125
126char *test_floorf(uint32 *in, uint32 *out) {
127 test_rintf(in, out, isfloor: 1, isceil: 0);
128 return NULL;
129}
130
131char *test_fmod(uint32 *a, uint32 *b, uint32 *out) {
132 int sign;
133 int32 aex, bex;
134 uint32 am[2], bm[2];
135
136 if (((a[0] & 0x7FFFFFFF) << 1) + !!a[1] > 0xFFE00000 ||
137 ((b[0] & 0x7FFFFFFF) << 1) + !!b[1] > 0xFFE00000) {
138 /* a or b is NaN: return QNaN, optionally with IVO */
139 uint32 an, bn;
140 out[0] = 0x7ff80000;
141 out[1] = 1;
142 an = ((a[0] & 0x7FFFFFFF) << 1) + !!a[1];
143 bn = ((b[0] & 0x7FFFFFFF) << 1) + !!b[1];
144 if ((an > 0xFFE00000 && an < 0xFFF00000) ||
145 (bn > 0xFFE00000 && bn < 0xFFF00000))
146 return "i"; /* at least one SNaN: IVO */
147 else
148 return NULL; /* no SNaNs, but at least 1 QNaN */
149 }
150 if ((b[0] & 0x7FFFFFFF) == 0 && b[1] == 0) { /* b==0: EDOM */
151 out[0] = 0x7ff80000;
152 out[1] = 1;
153 return "EDOM status=i";
154 }
155 if ((a[0] & 0x7FF00000) == 0x7FF00000) { /* a==Inf: EDOM */
156 out[0] = 0x7ff80000;
157 out[1] = 1;
158 return "EDOM status=i";
159 }
160 if ((b[0] & 0x7FF00000) == 0x7FF00000) { /* b==Inf: return a */
161 out[0] = a[0];
162 out[1] = a[1];
163 return NULL;
164 }
165 if ((a[0] & 0x7FFFFFFF) == 0 && a[1] == 0) { /* a==0: return a */
166 out[0] = a[0];
167 out[1] = a[1];
168 return NULL;
169 }
170
171 /*
172 * OK. That's the special cases cleared out of the way. Now we
173 * have finite (though not necessarily normal) a and b.
174 */
175 sign = a[0] & 0x80000000; /* we discard sign of b */
176 test_frexp(x: a, out: am, nout: (uint32 *)&aex);
177 test_frexp(x: b, out: bm, nout: (uint32 *)&bex);
178 am[0] &= 0xFFFFF, am[0] |= 0x100000;
179 bm[0] &= 0xFFFFF, bm[0] |= 0x100000;
180
181 while (aex >= bex) {
182 if (am[0] > bm[0] || (am[0] == bm[0] && am[1] >= bm[1])) {
183 am[1] -= bm[1];
184 am[0] = am[0] - bm[0] - (am[1] > ~bm[1]);
185 }
186 if (aex > bex) {
187 am[0] = (am[0] << 1) | ((am[1] & 0x80000000) >> 31);
188 am[1] <<= 1;
189 aex--;
190 } else
191 break;
192 }
193
194 /*
195 * Renormalise final result; this can be cunningly done by
196 * passing a denormal to ldexp.
197 */
198 aex += 0x3fd;
199 am[0] |= sign;
200 test_ldexp(x: am, n: (uint32 *)&aex, out);
201
202 return NULL; /* FIXME */
203}
204
205char *test_fmodf(uint32 *a, uint32 *b, uint32 *out) {
206 int sign;
207 int32 aex, bex;
208 uint32 am, bm;
209
210 if ((*a & 0x7FFFFFFF) > 0x7F800000 ||
211 (*b & 0x7FFFFFFF) > 0x7F800000) {
212 /* a or b is NaN: return QNaN, optionally with IVO */
213 uint32 an, bn;
214 *out = 0x7fc00001;
215 an = *a & 0x7FFFFFFF;
216 bn = *b & 0x7FFFFFFF;
217 if ((an > 0x7f800000 && an < 0x7fc00000) ||
218 (bn > 0x7f800000 && bn < 0x7fc00000))
219 return "i"; /* at least one SNaN: IVO */
220 else
221 return NULL; /* no SNaNs, but at least 1 QNaN */
222 }
223 if ((*b & 0x7FFFFFFF) == 0) { /* b==0: EDOM */
224 *out = 0x7fc00001;
225 return "EDOM status=i";
226 }
227 if ((*a & 0x7F800000) == 0x7F800000) { /* a==Inf: EDOM */
228 *out = 0x7fc00001;
229 return "EDOM status=i";
230 }
231 if ((*b & 0x7F800000) == 0x7F800000) { /* b==Inf: return a */
232 *out = *a;
233 return NULL;
234 }
235 if ((*a & 0x7FFFFFFF) == 0) { /* a==0: return a */
236 *out = *a;
237 return NULL;
238 }
239
240 /*
241 * OK. That's the special cases cleared out of the way. Now we
242 * have finite (though not necessarily normal) a and b.
243 */
244 sign = a[0] & 0x80000000; /* we discard sign of b */
245 test_frexpf(x: a, out: &am, nout: (uint32 *)&aex);
246 test_frexpf(x: b, out: &bm, nout: (uint32 *)&bex);
247 am &= 0x7FFFFF, am |= 0x800000;
248 bm &= 0x7FFFFF, bm |= 0x800000;
249
250 while (aex >= bex) {
251 if (am >= bm) {
252 am -= bm;
253 }
254 if (aex > bex) {
255 am <<= 1;
256 aex--;
257 } else
258 break;
259 }
260
261 /*
262 * Renormalise final result; this can be cunningly done by
263 * passing a denormal to ldexp.
264 */
265 aex += 0x7d;
266 am |= sign;
267 test_ldexpf(x: &am, n: (uint32 *)&aex, out);
268
269 return NULL; /* FIXME */
270}
271
272char *test_ldexp(uint32 *x, uint32 *np, uint32 *out) {
273 int n = *np;
274 int32 n2;
275 uint32 y[2];
276 int ex = (x[0] >> 20) & 0x7FF; /* exponent */
277 int sign = x[0] & 0x80000000;
278
279 if (ex == 0x7FF) { /* inf/NaN; just return x */
280 out[0] = x[0];
281 out[1] = x[1];
282 return NULL;
283 }
284 if ((x[0] & 0x7FFFFFFF) == 0 && x[1] == 0) { /* zero: return x */
285 out[0] = x[0];
286 out[1] = x[1];
287 return NULL;
288 }
289
290 test_frexp(x, out: y, nout: (uint32 *)&n2);
291 ex = n + n2;
292 if (ex > 0x400) { /* overflow */
293 out[0] = sign | 0x7FF00000;
294 out[1] = 0;
295 return "overflow";
296 }
297 /*
298 * Underflow. 2^-1074 is 00000000.00000001; so if ex == -1074
299 * then we have something [2^-1075,2^-1074). Under round-to-
300 * nearest-even, this whole interval rounds up to 2^-1074,
301 * except for the bottom endpoint which rounds to even and is
302 * an underflow condition.
303 *
304 * So, ex < -1074 is definite underflow, and ex == -1074 is
305 * underflow iff all mantissa bits are zero.
306 */
307 if (ex < -1074 || (ex == -1074 && (y[0] & 0xFFFFF) == 0 && y[1] == 0)) {
308 out[0] = sign; /* underflow: correctly signed zero */
309 out[1] = 0;
310 return "underflow";
311 }
312
313 /*
314 * No overflow or underflow; should be nice and simple, unless
315 * we have to denormalise and round the result.
316 */
317 if (ex < -1021) { /* denormalise and round */
318 uint32 roundword;
319 y[0] &= 0x000FFFFF;
320 y[0] |= 0x00100000; /* set leading bit */
321 roundword = 0;
322 while (ex < -1021) {
323 if (roundword & 1)
324 roundword |= 2; /* preserve sticky bit */
325 roundword = (roundword >> 1) | ((y[1] & 1) << 31);
326 y[1] = (y[1] >> 1) | ((y[0] & 1) << 31);
327 y[0] = y[0] >> 1;
328 ex++;
329 }
330 if (roundword > 0x80000000 || /* round up */
331 (roundword == 0x80000000 && (y[1] & 1))) { /* round up to even */
332 y[1]++;
333 y[0] += (y[1] == 0);
334 }
335 out[0] = sign | y[0];
336 out[1] = y[1];
337 /* Proper ERANGE underflow was handled earlier, but we still
338 * expect an IEEE Underflow exception if this partially
339 * underflowed result is not exact. */
340 if (roundword)
341 return "u";
342 return NULL; /* underflow was handled earlier */
343 } else {
344 out[0] = y[0] + (ex << 20);
345 out[1] = y[1];
346 return NULL;
347 }
348}
349
350char *test_ldexpf(uint32 *x, uint32 *np, uint32 *out) {
351 int n = *np;
352 int32 n2;
353 uint32 y;
354 int ex = (*x >> 23) & 0xFF; /* exponent */
355 int sign = *x & 0x80000000;
356
357 if (ex == 0xFF) { /* inf/NaN; just return x */
358 *out = *x;
359 return NULL;
360 }
361 if ((*x & 0x7FFFFFFF) == 0) { /* zero: return x */
362 *out = *x;
363 return NULL;
364 }
365
366 test_frexpf(x, out: &y, nout: (uint32 *)&n2);
367 ex = n + n2;
368 if (ex > 0x80) { /* overflow */
369 *out = sign | 0x7F800000;
370 return "overflow";
371 }
372 /*
373 * Underflow. 2^-149 is 00000001; so if ex == -149 then we have
374 * something [2^-150,2^-149). Under round-to- nearest-even,
375 * this whole interval rounds up to 2^-149, except for the
376 * bottom endpoint which rounds to even and is an underflow
377 * condition.
378 *
379 * So, ex < -149 is definite underflow, and ex == -149 is
380 * underflow iff all mantissa bits are zero.
381 */
382 if (ex < -149 || (ex == -149 && (y & 0x7FFFFF) == 0)) {
383 *out = sign; /* underflow: correctly signed zero */
384 return "underflow";
385 }
386
387 /*
388 * No overflow or underflow; should be nice and simple, unless
389 * we have to denormalise and round the result.
390 */
391 if (ex < -125) { /* denormalise and round */
392 uint32 roundword;
393 y &= 0x007FFFFF;
394 y |= 0x00800000; /* set leading bit */
395 roundword = 0;
396 while (ex < -125) {
397 if (roundword & 1)
398 roundword |= 2; /* preserve sticky bit */
399 roundword = (roundword >> 1) | ((y & 1) << 31);
400 y = y >> 1;
401 ex++;
402 }
403 if (roundword > 0x80000000 || /* round up */
404 (roundword == 0x80000000 && (y & 1))) { /* round up to even */
405 y++;
406 }
407 *out = sign | y;
408 /* Proper ERANGE underflow was handled earlier, but we still
409 * expect an IEEE Underflow exception if this partially
410 * underflowed result is not exact. */
411 if (roundword)
412 return "u";
413 return NULL; /* underflow was handled earlier */
414 } else {
415 *out = y + (ex << 23);
416 return NULL;
417 }
418}
419
420char *test_frexp(uint32 *x, uint32 *out, uint32 *nout) {
421 int ex = (x[0] >> 20) & 0x7FF; /* exponent */
422 if (ex == 0x7FF) { /* inf/NaN; return x/0 */
423 out[0] = x[0];
424 out[1] = x[1];
425 nout[0] = 0;
426 return NULL;
427 }
428 if (ex == 0) { /* denormals/zeros */
429 int sign;
430 uint32 xh, xl;
431 if ((x[0] & 0x7FFFFFFF) == 0 && x[1] == 0) {
432 /* zero: return x/0 */
433 out[0] = x[0];
434 out[1] = x[1];
435 nout[0] = 0;
436 return NULL;
437 }
438 sign = x[0] & 0x80000000;
439 xh = x[0] & 0x7FFFFFFF;
440 xl = x[1];
441 ex = 1;
442 while (!(xh & 0x100000)) {
443 ex--;
444 xh = (xh << 1) | ((xl >> 31) & 1);
445 xl = (xl & 0x7FFFFFFF) << 1;
446 }
447 out[0] = sign | 0x3FE00000 | (xh & 0xFFFFF);
448 out[1] = xl;
449 nout[0] = ex - 0x3FE;
450 return NULL;
451 }
452 out[0] = 0x3FE00000 | (x[0] & 0x800FFFFF);
453 out[1] = x[1];
454 nout[0] = ex - 0x3FE;
455 return NULL; /* ordinary number; no error */
456}
457
458char *test_frexpf(uint32 *x, uint32 *out, uint32 *nout) {
459 int ex = (*x >> 23) & 0xFF; /* exponent */
460 if (ex == 0xFF) { /* inf/NaN; return x/0 */
461 *out = *x;
462 nout[0] = 0;
463 return NULL;
464 }
465 if (ex == 0) { /* denormals/zeros */
466 int sign;
467 uint32 xv;
468 if ((*x & 0x7FFFFFFF) == 0) {
469 /* zero: return x/0 */
470 *out = *x;
471 nout[0] = 0;
472 return NULL;
473 }
474 sign = *x & 0x80000000;
475 xv = *x & 0x7FFFFFFF;
476 ex = 1;
477 while (!(xv & 0x800000)) {
478 ex--;
479 xv = xv << 1;
480 }
481 *out = sign | 0x3F000000 | (xv & 0x7FFFFF);
482 nout[0] = ex - 0x7E;
483 return NULL;
484 }
485 *out = 0x3F000000 | (*x & 0x807FFFFF);
486 nout[0] = ex - 0x7E;
487 return NULL; /* ordinary number; no error */
488}
489
490char *test_modf(uint32 *x, uint32 *fout, uint32 *iout) {
491 int ex = (x[0] >> 20) & 0x7FF; /* exponent */
492 int sign = x[0] & 0x80000000;
493 uint32 fh, fl;
494
495 if (((x[0] & 0x7FFFFFFF) | (!!x[1])) > 0x7FF00000) {
496 /*
497 * NaN input: return the same in _both_ outputs.
498 */
499 fout[0] = iout[0] = x[0];
500 fout[1] = iout[1] = x[1];
501 return NULL;
502 }
503
504 test_rint(in: x, out: iout, isfloor: 0, isceil: 0);
505 fh = x[0] - iout[0];
506 fl = x[1] - iout[1];
507 if (!fh && !fl) { /* no fraction part */
508 fout[0] = sign;
509 fout[1] = 0;
510 return NULL;
511 }
512 if (!(iout[0] & 0x7FFFFFFF) && !iout[1]) { /* no integer part */
513 fout[0] = x[0];
514 fout[1] = x[1];
515 return NULL;
516 }
517 while (!(fh & 0x100000)) {
518 ex--;
519 fh = (fh << 1) | ((fl >> 31) & 1);
520 fl = (fl & 0x7FFFFFFF) << 1;
521 }
522 fout[0] = sign | (ex << 20) | (fh & 0xFFFFF);
523 fout[1] = fl;
524 return NULL;
525}
526
527char *test_modff(uint32 *x, uint32 *fout, uint32 *iout) {
528 int ex = (*x >> 23) & 0xFF; /* exponent */
529 int sign = *x & 0x80000000;
530 uint32 f;
531
532 if ((*x & 0x7FFFFFFF) > 0x7F800000) {
533 /*
534 * NaN input: return the same in _both_ outputs.
535 */
536 *fout = *iout = *x;
537 return NULL;
538 }
539
540 test_rintf(in: x, out: iout, isfloor: 0, isceil: 0);
541 f = *x - *iout;
542 if (!f) { /* no fraction part */
543 *fout = sign;
544 return NULL;
545 }
546 if (!(*iout & 0x7FFFFFFF)) { /* no integer part */
547 *fout = *x;
548 return NULL;
549 }
550 while (!(f & 0x800000)) {
551 ex--;
552 f = f << 1;
553 }
554 *fout = sign | (ex << 23) | (f & 0x7FFFFF);
555 return NULL;
556}
557
558char *test_copysign(uint32 *x, uint32 *y, uint32 *out)
559{
560 int ysign = y[0] & 0x80000000;
561 int xhigh = x[0] & 0x7fffffff;
562
563 out[0] = ysign | xhigh;
564 out[1] = x[1];
565
566 /* There can be no error */
567 return NULL;
568}
569
570char *test_copysignf(uint32 *x, uint32 *y, uint32 *out)
571{
572 int ysign = y[0] & 0x80000000;
573 int xhigh = x[0] & 0x7fffffff;
574
575 out[0] = ysign | xhigh;
576
577 /* There can be no error */
578 return NULL;
579}
580
581char *test_isfinite(uint32 *x, uint32 *out)
582{
583 int xhigh = x[0];
584 /* Being finite means that the exponent is not 0x7ff */
585 if ((xhigh & 0x7ff00000) == 0x7ff00000) out[0] = 0;
586 else out[0] = 1;
587 return NULL;
588}
589
590char *test_isfinitef(uint32 *x, uint32 *out)
591{
592 /* Being finite means that the exponent is not 0xff */
593 if ((x[0] & 0x7f800000) == 0x7f800000) out[0] = 0;
594 else out[0] = 1;
595 return NULL;
596}
597
598char *test_isinff(uint32 *x, uint32 *out)
599{
600 /* Being infinite means that our bottom 30 bits equate to 0x7f800000 */
601 if ((x[0] & 0x7fffffff) == 0x7f800000) out[0] = 1;
602 else out[0] = 0;
603 return NULL;
604}
605
606char *test_isinf(uint32 *x, uint32 *out)
607{
608 int xhigh = x[0];
609 int xlow = x[1];
610 /* Being infinite means that our fraction is zero and exponent is 0x7ff */
611 if (((xhigh & 0x7fffffff) == 0x7ff00000) && (xlow == 0)) out[0] = 1;
612 else out[0] = 0;
613 return NULL;
614}
615
616char *test_isnanf(uint32 *x, uint32 *out)
617{
618 /* Being NaN means that our exponent is 0xff and non-0 fraction */
619 int exponent = x[0] & 0x7f800000;
620 int fraction = x[0] & 0x007fffff;
621 if ((exponent == 0x7f800000) && (fraction != 0)) out[0] = 1;
622 else out[0] = 0;
623 return NULL;
624}
625
626char *test_isnan(uint32 *x, uint32 *out)
627{
628 /* Being NaN means that our exponent is 0x7ff and non-0 fraction */
629 int exponent = x[0] & 0x7ff00000;
630 int fractionhigh = x[0] & 0x000fffff;
631 if ((exponent == 0x7ff00000) && ((fractionhigh != 0) || x[1] != 0))
632 out[0] = 1;
633 else out[0] = 0;
634 return NULL;
635}
636
637char *test_isnormalf(uint32 *x, uint32 *out)
638{
639 /* Being normal means exponent is not 0 and is not 0xff */
640 int exponent = x[0] & 0x7f800000;
641 if (exponent == 0x7f800000) out[0] = 0;
642 else if (exponent == 0) out[0] = 0;
643 else out[0] = 1;
644 return NULL;
645}
646
647char *test_isnormal(uint32 *x, uint32 *out)
648{
649 /* Being normal means exponent is not 0 and is not 0x7ff */
650 int exponent = x[0] & 0x7ff00000;
651 if (exponent == 0x7ff00000) out[0] = 0;
652 else if (exponent == 0) out[0] = 0;
653 else out[0] = 1;
654 return NULL;
655}
656
657char *test_signbitf(uint32 *x, uint32 *out)
658{
659 /* Sign bit is bit 31 */
660 out[0] = (x[0] >> 31) & 1;
661 return NULL;
662}
663
664char *test_signbit(uint32 *x, uint32 *out)
665{
666 /* Sign bit is bit 31 */
667 out[0] = (x[0] >> 31) & 1;
668 return NULL;
669}
670
671char *test_fpclassify(uint32 *x, uint32 *out)
672{
673 int exponent = (x[0] & 0x7ff00000) >> 20;
674 int fraction = (x[0] & 0x000fffff) | x[1];
675
676 if ((exponent == 0x00) && (fraction == 0)) out[0] = 0;
677 else if ((exponent == 0x00) && (fraction != 0)) out[0] = 4;
678 else if ((exponent == 0x7ff) && (fraction == 0)) out[0] = 3;
679 else if ((exponent == 0x7ff) && (fraction != 0)) out[0] = 7;
680 else out[0] = 5;
681 return NULL;
682}
683
684char *test_fpclassifyf(uint32 *x, uint32 *out)
685{
686 int exponent = (x[0] & 0x7f800000) >> 23;
687 int fraction = x[0] & 0x007fffff;
688
689 if ((exponent == 0x000) && (fraction == 0)) out[0] = 0;
690 else if ((exponent == 0x000) && (fraction != 0)) out[0] = 4;
691 else if ((exponent == 0xff) && (fraction == 0)) out[0] = 3;
692 else if ((exponent == 0xff) && (fraction != 0)) out[0] = 7;
693 else out[0] = 5;
694 return NULL;
695}
696
697/*
698 * Internal function that compares doubles in x & y and returns -3, -2, -1, 0,
699 * 1 if they compare to be signaling, unordered, less than, equal or greater
700 * than.
701 */
702static int fpcmp4(uint32 *x, uint32 *y)
703{
704 int result = 0;
705
706 /*
707 * Sort out whether results are ordered or not to begin with
708 * NaNs have exponent 0x7ff, and non-zero fraction. Signaling NaNs take
709 * higher priority than quiet ones.
710 */
711 if ((x[0] & 0x7fffffff) >= 0x7ff80000) result = -2;
712 else if ((x[0] & 0x7fffffff) > 0x7ff00000) result = -3;
713 else if (((x[0] & 0x7fffffff) == 0x7ff00000) && (x[1] != 0)) result = -3;
714 if ((y[0] & 0x7fffffff) >= 0x7ff80000 && result != -3) result = -2;
715 else if ((y[0] & 0x7fffffff) > 0x7ff00000) result = -3;
716 else if (((y[0] & 0x7fffffff) == 0x7ff00000) && (y[1] != 0)) result = -3;
717 if (result != 0) return result;
718
719 /*
720 * The two forms of zero are equal
721 */
722 if (((x[0] & 0x7fffffff) == 0) && x[1] == 0 &&
723 ((y[0] & 0x7fffffff) == 0) && y[1] == 0)
724 return 0;
725
726 /*
727 * If x and y have different signs we can tell that they're not equal
728 * If x is +ve we have x > y return 1 - otherwise y is +ve return -1
729 */
730 if ((x[0] >> 31) != (y[0] >> 31))
731 return ((x[0] >> 31) == 0) - ((y[0] >> 31) == 0);
732
733 /*
734 * Now we have both signs the same, let's do an initial compare of the
735 * values.
736 *
737 * Whoever designed IEEE754's floating point formats is very clever and
738 * earns my undying admiration. Once you remove the sign-bit, the
739 * floating point numbers can be ordered using the standard <, ==, >
740 * operators will treating the fp-numbers as integers with that bit-
741 * pattern.
742 */
743 if ((x[0] & 0x7fffffff) < (y[0] & 0x7fffffff)) result = -1;
744 else if ((x[0] & 0x7fffffff) > (y[0] & 0x7fffffff)) result = 1;
745 else if (x[1] < y[1]) result = -1;
746 else if (x[1] > y[1]) result = 1;
747 else result = 0;
748
749 /*
750 * Now we return the result - is x is positive (and therefore so is y) we
751 * return the plain result - otherwise we negate it and return.
752 */
753 if ((x[0] >> 31) == 0) return result;
754 else return -result;
755}
756
757/*
758 * Internal function that compares floats in x & y and returns -3, -2, -1, 0,
759 * 1 if they compare to be signaling, unordered, less than, equal or greater
760 * than.
761 */
762static int fpcmp4f(uint32 *x, uint32 *y)
763{
764 int result = 0;
765
766 /*
767 * Sort out whether results are ordered or not to begin with
768 * NaNs have exponent 0xff, and non-zero fraction - we have to handle all
769 * signaling cases over the quiet ones
770 */
771 if ((x[0] & 0x7fffffff) >= 0x7fc00000) result = -2;
772 else if ((x[0] & 0x7fffffff) > 0x7f800000) result = -3;
773 if ((y[0] & 0x7fffffff) >= 0x7fc00000 && result != -3) result = -2;
774 else if ((y[0] & 0x7fffffff) > 0x7f800000) result = -3;
775 if (result != 0) return result;
776
777 /*
778 * The two forms of zero are equal
779 */
780 if (((x[0] & 0x7fffffff) == 0) && ((y[0] & 0x7fffffff) == 0))
781 return 0;
782
783 /*
784 * If x and y have different signs we can tell that they're not equal
785 * If x is +ve we have x > y return 1 - otherwise y is +ve return -1
786 */
787 if ((x[0] >> 31) != (y[0] >> 31))
788 return ((x[0] >> 31) == 0) - ((y[0] >> 31) == 0);
789
790 /*
791 * Now we have both signs the same, let's do an initial compare of the
792 * values.
793 *
794 * Whoever designed IEEE754's floating point formats is very clever and
795 * earns my undying admiration. Once you remove the sign-bit, the
796 * floating point numbers can be ordered using the standard <, ==, >
797 * operators will treating the fp-numbers as integers with that bit-
798 * pattern.
799 */
800 if ((x[0] & 0x7fffffff) < (y[0] & 0x7fffffff)) result = -1;
801 else if ((x[0] & 0x7fffffff) > (y[0] & 0x7fffffff)) result = 1;
802 else result = 0;
803
804 /*
805 * Now we return the result - is x is positive (and therefore so is y) we
806 * return the plain result - otherwise we negate it and return.
807 */
808 if ((x[0] >> 31) == 0) return result;
809 else return -result;
810}
811
812char *test_isgreater(uint32 *x, uint32 *y, uint32 *out)
813{
814 int result = fpcmp4(x, y);
815 *out = (result == 1);
816 return result == -3 ? "i" : NULL;
817}
818
819char *test_isgreaterequal(uint32 *x, uint32 *y, uint32 *out)
820{
821 int result = fpcmp4(x, y);
822 *out = (result >= 0);
823 return result == -3 ? "i" : NULL;
824}
825
826char *test_isless(uint32 *x, uint32 *y, uint32 *out)
827{
828 int result = fpcmp4(x, y);
829 *out = (result == -1);
830 return result == -3 ? "i" : NULL;
831}
832
833char *test_islessequal(uint32 *x, uint32 *y, uint32 *out)
834{
835 int result = fpcmp4(x, y);
836 *out = (result == -1) || (result == 0);
837 return result == -3 ? "i" : NULL;
838}
839
840char *test_islessgreater(uint32 *x, uint32 *y, uint32 *out)
841{
842 int result = fpcmp4(x, y);
843 *out = (result == -1) || (result == 1);
844 return result == -3 ? "i" : NULL;
845}
846
847char *test_isunordered(uint32 *x, uint32 *y, uint32 *out)
848{
849 int normal = 0;
850 int result = fpcmp4(x, y);
851
852 test_isnormal(x, out);
853 normal |= *out;
854 test_isnormal(x: y, out);
855 normal |= *out;
856 *out = (result == -2) || (result == -3);
857 return result == -3 ? "i" : NULL;
858}
859
860char *test_isgreaterf(uint32 *x, uint32 *y, uint32 *out)
861{
862 int result = fpcmp4f(x, y);
863 *out = (result == 1);
864 return result == -3 ? "i" : NULL;
865}
866
867char *test_isgreaterequalf(uint32 *x, uint32 *y, uint32 *out)
868{
869 int result = fpcmp4f(x, y);
870 *out = (result >= 0);
871 return result == -3 ? "i" : NULL;
872}
873
874char *test_islessf(uint32 *x, uint32 *y, uint32 *out)
875{
876 int result = fpcmp4f(x, y);
877 *out = (result == -1);
878 return result == -3 ? "i" : NULL;
879}
880
881char *test_islessequalf(uint32 *x, uint32 *y, uint32 *out)
882{
883 int result = fpcmp4f(x, y);
884 *out = (result == -1) || (result == 0);
885 return result == -3 ? "i" : NULL;
886}
887
888char *test_islessgreaterf(uint32 *x, uint32 *y, uint32 *out)
889{
890 int result = fpcmp4f(x, y);
891 *out = (result == -1) || (result == 1);
892 return result == -3 ? "i" : NULL;
893}
894
895char *test_isunorderedf(uint32 *x, uint32 *y, uint32 *out)
896{
897 int normal = 0;
898 int result = fpcmp4f(x, y);
899
900 test_isnormalf(x, out);
901 normal |= *out;
902 test_isnormalf(x: y, out);
903 normal |= *out;
904 *out = (result == -2) || (result == -3);
905 return result == -3 ? "i" : NULL;
906}
907

source code of libc/AOR_v20.02/math/test/rtest/semi.c