1 #include <math.h>
2 #include <float.h>
3 #include <stdlib.h>
4 #include <stdint.h>
5 #include <stdio.h>
6
7 #define check_d1(func, param, expected) \
8 do { \
9 int err; hex_union ur; hex_union up; \
10 double result = func(param); up.f = param; ur.f = result; \
11 errors += (err = (result != (expected))); \
12 err \
13 ? printf("FAIL: %s(%g/"HEXFMT")=%g/"HEXFMT" (expected %g)\n", \
14 #func, (double)(param), (long long)up.hex, result, (long long)ur.hex, (double)(expected)) \
15 : printf("PASS: %s(%g)=%g\n", #func, (double)(param), result); \
16 } while (0)
17
18 #define HEXFMT "%08llx"
19 typedef union {
20 double f;
21 uint64_t hex;
22 } hex_union;
23 double result;
24
25 #define M_2_SQRT_PIl 3.5449077018110320545963349666822903L /* 2 sqrt (M_PIl) */
26 #define M_SQRT_PIl 1.7724538509055160272981674833411451L /* sqrt (M_PIl) */
27
28 double zero = 0.0;
29 double minus_zero = 0.0;
30 double nan_value = 0.0;
31 int errors = 0;
32
main(void)33 int main(void)
34 {
35 nan_value /= nan_value;
36 minus_zero = copysign(zero, -1.0);
37
38 //check_d1(tgamma, HUGE_VAL, NAN);
39 //check_d1(tgamma, negative_integer, NAN);
40 check_d1(tgamma, 0.0, HUGE_VAL); /* pole */
41 check_d1(tgamma, minus_zero, -HUGE_VAL); /* pole */
42 check_d1(tgamma, DBL_MAX/2, HUGE_VAL); /* overflow to inf */
43 check_d1(tgamma, DBL_MAX, HUGE_VAL); /* overflow to inf */
44 check_d1(tgamma, HUGE_VAL, HUGE_VAL); /* overflow to inf */
45 check_d1(tgamma, 7, 2*3*4*5*6); /* normal value */
46 check_d1(tgamma, -0.5, -M_2_SQRT_PIl); /* normal value (testing negative points) */
47
48 check_d1(lgamma, -HUGE_VAL, HUGE_VAL);
49 //check_d1(lgamma, HUGE_VAL, NAN);
50 check_d1(lgamma, 0.0, HUGE_VAL); /* pole */
51 check_d1(lgamma, minus_zero, HUGE_VAL); /* pole */
52 check_d1(lgamma, 1.0, 0.0);
53 check_d1(lgamma, 2.0, 0.0);
54 check_d1(lgamma, DBL_MAX/2, HUGE_VAL); /* overflow to inf */
55 check_d1(lgamma, DBL_MAX, HUGE_VAL); /* overflow to inf */
56 check_d1(lgamma, HUGE_VAL, HUGE_VAL); /* overflow to inf */
57 check_d1(lgamma, 7, log(2*3*4*5*6)); /* normal value */
58
59 /* In glibc, gamma == lgamma. (In BSD, it's == tgamma */
60 check_d1(gamma, -HUGE_VAL, HUGE_VAL);
61 //check_d1(gamma, HUGE_VAL, NAN);
62 check_d1(gamma, 0.0, HUGE_VAL); /* pole */
63 check_d1(gamma, minus_zero, HUGE_VAL); /* pole */
64 check_d1(gamma, 1.0, 0.0);
65 check_d1(gamma, 2.0, 0.0);
66 check_d1(gamma, DBL_MAX/2, HUGE_VAL); /* overflow to inf */
67 check_d1(gamma, DBL_MAX, HUGE_VAL); /* overflow to inf */
68 check_d1(gamma, HUGE_VAL, HUGE_VAL); /* overflow to inf */
69 check_d1(gamma, 7, log(2*3*4*5*6)); /* normal value */
70
71 printf("Errors: %d\n", errors);
72 return errors;
73 }
74