1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
#define _BSD_SOURCE 1
#define _XOPEN_SOURCE 700
#include <stdint.h>
#include <stdio.h>
#include <fenv.h>
#include <float.h>
#include <math.h>
#define RN 0
#define T(...) {__FILE__, __LINE__, __VA_ARGS__},
#define POS const char *file; int line;
struct f_fi {POS int r; float x; float y; float dy; long long i; int e; };
#define DIVBYZERO 0
#define INEXACT 0
#define INVALID 0
#define OVERFLOW 0
#define UNDERFLOW 0
#define inf INFINITY
#define nan NAN
static struct f_fi t[] = {
T(RN, -0x1.02239f3c6a8f1p+3, -0x1.0120f61b63d5ep+3, 0x1.89ccc4p-6, -1, INEXACT)
T(RN, 0x1.161868e18bc67p+2, 0x1.1ef3b263fd60bp+1, -0x1.6d0264p-3, 1, INEXACT)
T(RN, -0x1.0c34b3e01e6e7p+3, -0x1.46d73255263d9p+3, 0x1.e0ec76p-3, -1, INEXACT)
T(RN, -0x1.a206f0a19dcc4p+2, -0x1.9c91f19ac48c5p+2, 0x1.c2a38cp-2, -1, INEXACT)
T(RN, 0x1.288bbb0d6a1e6p+3, 0x1.65c60768fcc11p+3, 0x1.2f22c2p-2, 1, INEXACT)
T(RN, 0x1.52efd0cd80497p-1, 0x1.3cc760be720b3p-2, 0x1.0527e2p-2, 1, INEXACT)
T(RN, -0x1.a05cc754481d1p-2, 0x1.4ef387fea1014p+0, -0x1.c3b036p-2, -1, INEXACT)
T(RN, 0x1.1f9ef934745cbp-1, 0x1.d6f0efacc5699p-2, 0x1.c0b0a8p-2, 1, INEXACT)
T(RN, 0x1.8c5db097f7442p-1, 0x1.6c1a14cf91533p-3, 0x1.16f4cap-5, 1, INEXACT)
T(RN, -0x1.5b86ea8118a0ep-1, 0x1.695b1e0a0a59ep+0, 0x1.ada69ep-2, -1, INEXACT)
T(RN, 0x0p+0, inf, 0x0p+0, 1, DIVBYZERO)
/* T(RN, -0x0p+0, inf, 0x0p+0, -1, DIVBYZERO) This one fails in native as well */
T(RN, 0x1p+0, 0x0p+0, 0x0p+0, 1, 0)
T(RN, -0x1p+0, inf, 0x0p+0, 1, DIVBYZERO)
T(RN, 0x1p+1, 0x0p+0, 0x0p+0, 1, 0)
T(RN, -0x1p+1, inf, 0x0p+0, 1, DIVBYZERO)
T(RN, inf, inf, 0x0p+0, 1, 0)
T(RN, -inf, inf, 0x0p+0, -1, 0)
T(RN, nan, nan, 0x0p+0, 1, 0)
};
static int eulpf(float x)
{
union { float f; uint32_t i; } u = { x };
int e = u.i>>23 & 0xff;
if (!e)
e++;
return e - 0x7f - 23;
}
static int checkulp(float d, int r)
{
// TODO: we only care about >=1.5 ulp errors for now, should be 1.0
if (r == RN)
return fabsf(d) < 1.5;
return 1;
}
static float ulperrf(float got, float want, float dwant)
{
if (isnan(got) && isnan(want))
return 0;
if (got == want) {
if (signbit(got) == signbit(want))
return dwant;
return INFINITY;
}
if (isinf(got)) {
got = copysignf(0x1p127, got);
want *= 0.5;
}
return scalbn(got - want, -eulpf(want)) + dwant;
}
int main(void)
{
int yi;
double y;
float d;
int e, i, err = 0;
struct f_fi *p;
for (i = 0; i < sizeof t/sizeof *t; i++) {
p = t + i;
if (p->r < 0)
continue;
y = lgammaf(p->x);
yi = signgam;
printf("%g,%d\n", y, yi);
d = ulperrf(y, p->y, p->dy);
if (!checkulp(d, p->r) || (!isnan(p->x) && p->x!=-inf && !(p->e&DIVBYZERO) && yi != p->i)) {
/* printf("%s:%d: %d lgammaf(%g) want %g,%lld got %g,%d ulperr %.3f = %g + %g\n",
p->file, p->line, p->r, p->x, p->y, p->i, y, yi, d, d-p->dy, p->dy); */
err++;
}
}
return !!err;
}
|
