non-nearest rounding ulp check
[libc-test] / src / common / mtest.c
1 #include <stdio.h>
2 #include <stdint.h>
3 #include "mtest.h"
4
5 int eulpf(float x)
6 {
7         union { float f; uint32_t i; } u = { x };
8         int e = u.i>>23 & 0xff;
9
10         if (!e)
11                 e++;
12         return e - 0x7f - 23;
13 }
14
15 int eulp(double x)
16 {
17         union { double f; uint64_t i; } u = { x };
18         int e = u.i>>52 & 0x7ff;
19
20         if (!e)
21                 e++;
22         return e - 0x3ff - 52;
23 }
24
25 int eulpl(long double x)
26 {
27 #if LDBL_MANT_DIG == 53
28         return eulp(x);
29 #elif LDBL_MANT_DIG == 64
30         union { long double f; struct {uint64_t m; uint16_t e; uint16_t pad;} i; } u = { x };
31         int e = u.i.e & 0x7fff;
32
33         if (!e)
34                 e++;
35         return e - 0x3fff - 63;
36 #else
37         // TODO
38         return 0;
39 #endif
40 }
41
42 float ulperrf(float got, float want, float dwant)
43 {
44         if (isnan(got) && isnan(want))
45                 return 0;
46         if (got == want) {
47                 if (signbit(got) == signbit(want))
48                         return dwant;
49                 return inf;
50         }
51         if (isinf(got)) {
52                 got = copysignf(0x1p127, got);
53                 want *= 0.5;
54         }
55         return scalbn(got - want, -eulpf(want)) + dwant;
56 }
57
58 float ulperr(double got, double want, float dwant)
59 {
60         if (isnan(got) && isnan(want))
61                 return 0;
62         if (got == want) {
63                 if (signbit(got) == signbit(want))
64                         return dwant;
65                 return inf; // treat 0 sign errors badly
66         }
67         if (isinf(got)) {
68                 got = copysign(0x1p1023, got);
69                 want *= 0.5;
70         }
71         return scalbn(got - want, -eulp(want)) + dwant;
72 }
73
74 float ulperrl(long double got, long double want, float dwant)
75 {
76 #if LDBL_MANT_DIG == 53
77         return ulperr(got, want, dwant);
78 #elif LDBL_MANT_DIG == 64
79         if (isnan(got) && isnan(want))
80                 return 0;
81         if (got == want) {
82                 if (signbit(got) == signbit(want))
83                         return dwant;
84                 return inf;
85         }
86         if (isinf(got)) {
87                 got = copysignl(0x1p16383L, got);
88                 want *= 0.5;
89         }
90         return scalbnl(got - want, -eulpl(want)) + dwant;
91 #else
92         // TODO
93         return inf;
94 #endif
95 }
96
97 #define length(a) (sizeof(a)/sizeof*(a))
98 #define flag(x) {x, #x}
99 static struct {
100         int flag;
101         char *s;
102 } eflags[] = {
103         flag(INEXACT),
104         flag(INVALID),
105         flag(DIVBYZERO),
106         flag(UNDERFLOW),
107         flag(OVERFLOW)
108 };
109
110 char *estr(int f)
111 {
112         static char buf[256];
113         char *p = buf;
114         int i, all = 0;
115
116         for (i = 0; i < length(eflags); i++)
117                 if (f & eflags[i].flag) {
118                         p += sprintf(p, "%s%s", all ? "|" : "", eflags[i].s);
119                         all |= eflags[i].flag;
120                 }
121         if (all != f) {
122                 p += sprintf(p, "%s%d", all ? "|" : "", f & ~all);
123                 all = f;
124         }
125         p += sprintf(p, "%s", all ? "" : "0");
126         return buf;
127 }
128
129 char *rstr(int r)
130 {
131         switch (r) {
132         case RN: return "RN";
133 #ifdef FE_TOWARDZERO
134         case RZ: return "RZ";
135 #endif
136 #ifdef FE_UPWARD
137         case RU: return "RU";
138 #endif
139 #ifdef FE_DOWNWARD
140         case RD: return "RD";
141 #endif
142         }
143         return "R?";
144 }