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