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