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
80static inline long double complex
81CMPLXL(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/*
92 * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0.
93 * Fail an assertion if they differ.
94 */
95static int
96fpequal(long double d1, long double d2)
97{
98
99	if (d1 != d2)
100		return (isnan(d1) && isnan(d2));
101	return (copysignl(1.0, d1) == copysignl(1.0, d2));
102}
103
104/*
105 * Determine whether x and y are equal, with two special rules:
106 *	+0.0 != -0.0
107 *	 NaN == NaN
108 * If checksign is 0, we compare the absolute values instead.
109 */
110static int
111fpequal_cs(long double x, long double y, int checksign)
112{
113	if (isnan(x) && isnan(y))
114		return (1);
115	if (checksign)
116		return (x == y && !signbit(x) == !signbit(y));
117	else
118		return (fabsl(x) == fabsl(y));
119}
120
121static int
122fpequal_tol(long double x, long double y, long double tol, unsigned int flags)
123{
124	fenv_t env;
125	int ret;
126
127	if (isnan(x) && isnan(y))
128		return (1);
129	if (!signbit(x) != !signbit(y) && (flags & CS_BOTH))
130		return (0);
131	if (x == y)
132		return (1);
133	if (tol == 0)
134		return (0);
135
136	/* Hard case: need to check the tolerance. */
137	feholdexcept(&env);
138	/*
139	 * For our purposes here, if y=0, we interpret tol as an absolute
140	 * tolerance. This is to account for roundoff in the input, e.g.,
141	 * cos(Pi/2) ~= 0.
142	 */
143	if ((flags & FPE_ABS_ZERO) && y == 0.0)
144		ret = fabsl(x - y) <= fabsl(tol);
145	else
146		ret = fabsl(x - y) <= fabsl(y * tol);
147	fesetenv(&env);
148	return (ret);
149}
150
151static int
152cfpequal(long double complex d1, long double complex d2)
153{
154
155	return (fpequal(creall(d1), creall(d2)) &&
156		fpequal(cimagl(d1), cimagl(d2)));
157}
158
159static int
160cfpequal_cs(long double complex x, long double complex y, int checksign)
161{
162	return (fpequal_cs(creal(x), creal(y), checksign)
163		&& fpequal_cs(cimag(x), cimag(y), checksign));
164}
165
166static int
167cfpequal_tol(long double complex x, long double complex y, long double tol,
168		     unsigned int flags)
169{
170	return (fpequal_tol(creal(x), creal(y), tol, flags)
171		&& fpequal_tol(cimag(x), cimag(y), tol, flags));
172}
173
174#endif /* _TEST_UTILS_H_ */
175