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 #else
37 // TODO
38 return 0;
39 #endif
40 }
41
ulperrf(float got,float want,float dwant)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
ulperr(double got,double want,float dwant)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
ulperrl(long double got,long double want,float dwant)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
estr(int f)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
rstr(int r)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 }
145