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