1 /*- 2 * Copyright (c) 2005-2013 David Schultz <das@FreeBSD.org> 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 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 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 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 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 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 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 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