1 #include <stdio.h>
2 #include <stdint.h>
3 #include "mtest.h"
4
eulpf(float x)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
eulp(double x)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
eulpl(long double x)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 #elif LDBL_MANT_DIG == 113
37 union { long double f; struct {uint64_t lo; uint32_t mid; uint16_t top; uint16_t se;} i; } u = { x };
38 int e = u.i.se & 0x7fff;
39
40 if (!e)
41 e++;
42 return e - 0x3fff - 112;
43
44 #else
45 // TODO
46 return 0;
47 #endif
48 }
49
ulperrf(float got,float want,float dwant)50 float ulperrf(float got, float want, float dwant)
51 {
52 if (isnan(got) && isnan(want))
53 return 0;
54 if (got == want) {
55 if (signbit(got) == signbit(want))
56 return dwant;
57 return inf;
58 }
59 if (isinf(got)) {
60 got = copysignf(0x1p127, got);
61 want *= 0.5;
62 }
63 return scalbn(got - want, -eulpf(want)) + dwant;
64 }
65
ulperr(double got,double want,float dwant)66 float ulperr(double got, double want, float dwant)
67 {
68 if (isnan(got) && isnan(want))
69 return 0;
70 if (got == want) {
71 if (signbit(got) == signbit(want))
72 return dwant;
73 return inf; // treat 0 sign errors badly
74 }
75 if (isinf(got)) {
76 got = copysign(0x1p1023, got);
77 want *= 0.5;
78 }
79 return scalbn(got - want, -eulp(want)) + dwant;
80 }
81
ulperrl(long double got,long double want,float dwant)82 float ulperrl(long double got, long double want, float dwant)
83 {
84 #if LDBL_MANT_DIG == 53
85 return ulperr(got, want, dwant);
86 #elif LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113
87 if (isnan(got) && isnan(want))
88 return 0;
89 if (got == want) {
90 if (signbit(got) == signbit(want))
91 return dwant;
92 return inf;
93 }
94 if (isinf(got)) {
95 got = copysignl(0x1p16383L, got);
96 want *= 0.5;
97 }
98 return scalbnl(got - want, -eulpl(want)) + dwant;
99 #else
100 // TODO
101 return inf;
102 #endif
103 }
104
105 #define length(a) (sizeof(a)/sizeof*(a))
106 #define flag(x) {x, #x}
107 static struct {
108 int flag;
109 char *s;
110 } eflags[] = {
111 flag(INEXACT),
112 flag(INVALID),
113 flag(DIVBYZERO),
114 flag(UNDERFLOW),
115 flag(OVERFLOW)
116 };
117
estr(int f)118 char *estr(int f)
119 {
120 static char buf[256];
121 char *p = buf;
122 int i, all = 0;
123
124 for (i = 0; i < length(eflags); i++)
125 if (f & eflags[i].flag) {
126 p += sprintf(p, "%s%s", all ? "|" : "", eflags[i].s);
127 all |= eflags[i].flag;
128 }
129 if (all != f) {
130 p += sprintf(p, "%s%d", all ? "|" : "", f & ~all);
131 all = f;
132 }
133 p += sprintf(p, "%s", all ? "" : "0");
134 return buf;
135 }
136
rstr(int r)137 char *rstr(int r)
138 {
139 switch (r) {
140 case RN: return "RN";
141 #ifdef FE_TOWARDZERO
142 case RZ: return "RZ";
143 #endif
144 #ifdef FE_UPWARD
145 case RU: return "RU";
146 #endif
147 #ifdef FE_DOWNWARD
148 case RD: return "RD";
149 #endif
150 }
151 return "R?";
152 }
153