1 /*-
2 * Copyright (c) 2005-2013 David Schultz <[email protected]>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 *
26 * $FreeBSD$
27 */
28
29 #ifndef _TEST_UTILS_H_
30 #define _TEST_UTILS_H_
31
32 #include <complex.h>
33 #include <fenv.h>
34
35 /*
36 * Implementations are permitted to define additional exception flags
37 * not specified in the standard, so it is not necessarily true that
38 * FE_ALL_EXCEPT == ALL_STD_EXCEPT.
39 */
40 #define ALL_STD_EXCEPT (FE_DIVBYZERO | FE_INEXACT | FE_INVALID | \
41 FE_OVERFLOW | FE_UNDERFLOW)
42 #define OPT_INVALID (ALL_STD_EXCEPT & ~FE_INVALID)
43 #define OPT_INEXACT (ALL_STD_EXCEPT & ~FE_INEXACT)
44 #define FLT_ULP() ldexpl(1.0, 1 - FLT_MANT_DIG)
45 #define DBL_ULP() ldexpl(1.0, 1 - DBL_MANT_DIG)
46 #define LDBL_ULP() ldexpl(1.0, 1 - LDBL_MANT_DIG)
47
48 /*
49 * Flags that control the behavior of various fpequal* functions.
50 * XXX This is messy due to merging various notions of "close enough"
51 * that are best suited for different functions.
52 *
53 * CS_REAL
54 * CS_IMAG
55 * CS_BOTH
56 * (cfpequal_cs, fpequal_tol, cfpequal_tol) Whether to check the sign of
57 * the real part of the result, the imaginary part, or both.
58 *
59 * FPE_ABS_ZERO
60 * (fpequal_tol, cfpequal_tol) If set, treats the tolerance as an absolute
61 * tolerance when the expected value is 0. This is useful when there is
62 * round-off error in the input, e.g., cos(Pi/2) ~= 0.
63 */
64 #define CS_REAL 0x01
65 #define CS_IMAG 0x02
66 #define CS_BOTH (CS_REAL | CS_IMAG)
67 #define FPE_ABS_ZERO 0x04
68
69 #ifdef DEBUG
70 #define debug(...) printf(__VA_ARGS__)
71 #else
72 #define debug(...) (void)0
73 #endif
74
75 /*
76 * XXX The ancient version of gcc in the base system doesn't support CMPLXL,
77 * but we can fake it most of the time.
78 */
79 #ifndef CMPLXL
80 static inline long double complex
CMPLXL(long double x,long double y)81 CMPLXL(long double x, long double y)
82 {
83 long double complex z;
84
85 __real__ z = x;
86 __imag__ z = y;
87 return (z);
88 }
89 #endif
90
91 static int fpequal(long double, long double) __used;
92 static int cfpequal(long double complex, long double complex) __used;
93 static int cfpequal_cs(long double complex, long double complex,
94 int) __used;
95 static int cfpequal_tol(long double complex, long double complex,
96 long double, unsigned int) __used;
97
98 /*
99 * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0.
100 * Fail an assertion if they differ.
101 */
102 static int
fpequal(long double d1,long double d2)103 fpequal(long double d1, long double d2)
104 {
105
106 if (d1 != d2)
107 return (isnan(d1) && isnan(d2));
108 return (copysignl(1.0, d1) == copysignl(1.0, d2));
109 }
110
111 /*
112 * Determine whether x and y are equal, with two special rules:
113 * +0.0 != -0.0
114 * NaN == NaN
115 * If checksign is 0, we compare the absolute values instead.
116 */
117 static int
fpequal_cs(long double x,long double y,int checksign)118 fpequal_cs(long double x, long double y, int checksign)
119 {
120 if (isnan(x) && isnan(y))
121 return (1);
122 if (checksign)
123 return (x == y && !signbit(x) == !signbit(y));
124 else
125 return (fabsl(x) == fabsl(y));
126 }
127
128 static int
fpequal_tol(long double x,long double y,long double tol,unsigned int flags)129 fpequal_tol(long double x, long double y, long double tol,
130 unsigned int flags)
131 {
132 fenv_t env;
133 int ret;
134
135 if (isnan(x) && isnan(y))
136 return (1);
137 if (!signbit(x) != !signbit(y) && (flags & CS_BOTH))
138 return (0);
139 if (x == y)
140 return (1);
141 if (tol == 0)
142 return (0);
143
144 /* Hard case: need to check the tolerance. */
145 feholdexcept(&env);
146 /*
147 * For our purposes here, if y=0, we interpret tol as an absolute
148 * tolerance. This is to account for roundoff in the input, e.g.,
149 * cos(Pi/2) ~= 0.
150 */
151 if ((flags & FPE_ABS_ZERO) && y == 0.0)
152 ret = fabsl(x - y) <= fabsl(tol);
153 else
154 ret = fabsl(x - y) <= fabsl(y * tol);
155 fesetenv(&env);
156 return (ret);
157 }
158
159 static int
cfpequal(long double complex d1,long double complex d2)160 cfpequal(long double complex d1, long double complex d2)
161 {
162
163 return (fpequal(creall(d1), creall(d2)) &&
164 fpequal(cimagl(d1), cimagl(d2)));
165 }
166
167 static int
cfpequal_cs(long double complex x,long double complex y,int checksign)168 cfpequal_cs(long double complex x, long double complex y, int checksign)
169 {
170 return (fpequal_cs(creal(x), creal(y), checksign)
171 && fpequal_cs(cimag(x), cimag(y), checksign));
172 }
173
174 static int
cfpequal_tol(long double complex x,long double complex y,long double tol,unsigned int flags)175 cfpequal_tol(long double complex x, long double complex y, long double tol,
176 unsigned int flags)
177 {
178 return (fpequal_tol(creal(x), creal(y), tol, flags)
179 && fpequal_tol(cimag(x), cimag(y), tol, flags));
180 }
181
182 #endif /* _TEST_UTILS_H_ */
183