diff options
Diffstat (limited to 'tests/math/lgamma.in')
| -rw-r--r-- | tests/math/lgamma.in | 105 |
1 files changed, 105 insertions, 0 deletions
diff --git a/tests/math/lgamma.in b/tests/math/lgamma.in new file mode 100644 index 00000000..e96f5610 --- /dev/null +++ b/tests/math/lgamma.in @@ -0,0 +1,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; +} |