xref: /freebsd/lib/msun/tests/trig_test.c (revision 1d386b48)
1 /*-
2  * Copyright (c) 2008 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 
27 /*
28  * Tests for corner cases in trigonometric functions. Some accuracy tests
29  * are included as well, but these are very basic sanity checks, not
30  * intended to be comprehensive.
31  *
32  * The program for generating representable numbers near multiples of pi is
33  * available at http://www.cs.berkeley.edu/~wkahan/testpi/ .
34  */
35 
36 #include <sys/cdefs.h>
37 #include <sys/param.h>
38 
39 #include <fenv.h>
40 #include <float.h>
41 #include <math.h>
42 #include <stdio.h>
43 
44 #include "test-utils.h"
45 
46 #pragma STDC FENV_ACCESS ON
47 
48 /*
49  * Test that a function returns the correct value and sets the
50  * exception flags correctly. The exceptmask specifies which
51  * exceptions we should check. We need to be lenient for several
52  * reasons, but mainly because on some architectures it's impossible
53  * to raise FE_OVERFLOW without raising FE_INEXACT.
54  *
55  * These are macros instead of functions so that assert provides more
56  * meaningful error messages.
57  *
58  * XXX The volatile here is to avoid gcc's bogus constant folding and work
59  *     around the lack of support for the FENV_ACCESS pragma.
60  */
61 #define	test(func, x, result, exceptmask, excepts)	do {		\
62 	volatile long double _d = x;					\
63 	ATF_CHECK(feclearexcept(FE_ALL_EXCEPT) == 0);			\
64 	CHECK_FPEQUAL((func)(_d), (result));			\
65 	CHECK_FP_EXCEPTIONS_MSG(excepts, exceptmask, "for %s(%s)",	\
66 	    #func, #x);							\
67 } while (0)
68 
69 #define	testall(prefix, x, result, exceptmask, excepts)	do {		\
70 	test(prefix, x, (double)result, exceptmask, excepts);		\
71 	test(prefix##f, x, (float)result, exceptmask, excepts);		\
72 	test(prefix##l, x, result, exceptmask, excepts);		\
73 } while (0)
74 
75 #define	testdf(prefix, x, result, exceptmask, excepts)	do {		\
76 	test(prefix, x, (double)result, exceptmask, excepts);		\
77 	test(prefix##f, x, (float)result, exceptmask, excepts);		\
78 } while (0)
79 
80 ATF_TC(special);
81 ATF_TC_HEAD(special, tc)
82 {
83 
84 	atf_tc_set_md_var(tc, "descr",
85  	    "test special cases in sin(), cos(), and tan()");
86 }
87 ATF_TC_BODY(special, tc)
88 {
89 
90 	/* Values at 0 should be exact. */
91 	testall(tan, 0.0, 0.0, ALL_STD_EXCEPT, 0);
92 	testall(tan, -0.0, -0.0, ALL_STD_EXCEPT, 0);
93 	testall(cos, 0.0, 1.0, ALL_STD_EXCEPT, 0);
94 	testall(cos, -0.0, 1.0, ALL_STD_EXCEPT, 0);
95 	testall(sin, 0.0, 0.0, ALL_STD_EXCEPT, 0);
96 	testall(sin, -0.0, -0.0, ALL_STD_EXCEPT, 0);
97 
98 	/* func(+-Inf) == NaN */
99 	testall(tan, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
100 	testall(sin, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
101 	testall(cos, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
102 	testall(tan, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
103 	testall(sin, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
104 	testall(cos, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
105 
106 	/* func(NaN) == NaN */
107 	testall(tan, NAN, NAN, ALL_STD_EXCEPT, 0);
108 	testall(sin, NAN, NAN, ALL_STD_EXCEPT, 0);
109 	testall(cos, NAN, NAN, ALL_STD_EXCEPT, 0);
110 }
111 
112 #ifndef __i386__
113 ATF_TC(reduction);
114 ATF_TC_HEAD(reduction, tc)
115 {
116 
117 	atf_tc_set_md_var(tc, "descr",
118  	    "tests to ensure argument reduction for large arguments is accurate");
119 }
120 ATF_TC_BODY(reduction, tc)
121 {
122 	/* floats very close to odd multiples of pi */
123 	static const float f_pi_odd[] = {
124 		85563208.0f,
125 		43998769152.0f,
126 		9.2763667655669323e+25f,
127 		1.5458357838905804e+29f,
128 	};
129 	/* doubles very close to odd multiples of pi */
130 	static const double d_pi_odd[] = {
131 		3.1415926535897931,
132 		91.106186954104004,
133 		642615.9188844458,
134 		3397346.5699258847,
135 		6134899525417045.0,
136 		3.0213551960457761e+43,
137 		1.2646209897993783e+295,
138 		6.2083625380677099e+307,
139 	};
140 	/* long doubles very close to odd multiples of pi */
141 #if LDBL_MANT_DIG == 64
142 	static const long double ld_pi_odd[] = {
143 		1.1891886960373841596e+101L,
144 		1.07999475322710967206e+2087L,
145 		6.522151627890431836e+2147L,
146 		8.9368974898260328229e+2484L,
147 		9.2961044110572205863e+2555L,
148 		4.90208421886578286e+3189L,
149 		1.5275546401232615884e+3317L,
150 		1.7227465626338900093e+3565L,
151 		2.4160090594000745334e+3808L,
152 		9.8477555741888350649e+4314L,
153 		1.6061597222105160737e+4326L,
154 	};
155 #endif
156 
157 	unsigned i;
158 
159 #if defined(__amd64__) && defined(__clang__) && __clang_major__ >= 7 && \
160     __clang_major__ < 10 && __FreeBSD_cc_version < 1300002
161 	atf_tc_expect_fail("test fails with clang 7-9 - bug 234040");
162 #endif
163 
164 	for (i = 0; i < nitems(f_pi_odd); i++) {
165 		ATF_CHECK(fabs(sinf(f_pi_odd[i])) < FLT_EPSILON);
166 		ATF_CHECK(cosf(f_pi_odd[i]) == -1.0);
167 		ATF_CHECK(fabs(tan(f_pi_odd[i])) < FLT_EPSILON);
168 
169 		ATF_CHECK(fabs(sinf(-f_pi_odd[i])) < FLT_EPSILON);
170 		ATF_CHECK(cosf(-f_pi_odd[i]) == -1.0);
171 		ATF_CHECK(fabs(tanf(-f_pi_odd[i])) < FLT_EPSILON);
172 
173 		ATF_CHECK(fabs(sinf(f_pi_odd[i] * 2)) < FLT_EPSILON);
174 		ATF_CHECK(cosf(f_pi_odd[i] * 2) == 1.0);
175 		ATF_CHECK(fabs(tanf(f_pi_odd[i] * 2)) < FLT_EPSILON);
176 
177 		ATF_CHECK(fabs(sinf(-f_pi_odd[i] * 2)) < FLT_EPSILON);
178 		ATF_CHECK(cosf(-f_pi_odd[i] * 2) == 1.0);
179 		ATF_CHECK(fabs(tanf(-f_pi_odd[i] * 2)) < FLT_EPSILON);
180 	}
181 
182 	for (i = 0; i < nitems(d_pi_odd); i++) {
183 		ATF_CHECK(fabs(sin(d_pi_odd[i])) < 2 * DBL_EPSILON);
184 		ATF_CHECK(cos(d_pi_odd[i]) == -1.0);
185 		ATF_CHECK(fabs(tan(d_pi_odd[i])) < 2 * DBL_EPSILON);
186 
187 		ATF_CHECK(fabs(sin(-d_pi_odd[i])) < 2 * DBL_EPSILON);
188 		ATF_CHECK(cos(-d_pi_odd[i]) == -1.0);
189 		ATF_CHECK(fabs(tan(-d_pi_odd[i])) < 2 * DBL_EPSILON);
190 
191 		ATF_CHECK(fabs(sin(d_pi_odd[i] * 2)) < 2 * DBL_EPSILON);
192 		ATF_CHECK(cos(d_pi_odd[i] * 2) == 1.0);
193 		ATF_CHECK(fabs(tan(d_pi_odd[i] * 2)) < 2 * DBL_EPSILON);
194 
195 		ATF_CHECK(fabs(sin(-d_pi_odd[i] * 2)) < 2 * DBL_EPSILON);
196 		ATF_CHECK(cos(-d_pi_odd[i] * 2) == 1.0);
197 		ATF_CHECK(fabs(tan(-d_pi_odd[i] * 2)) < 2 * DBL_EPSILON);
198 	}
199 
200 #if LDBL_MANT_DIG == 64 /* XXX: || LDBL_MANT_DIG == 113 */
201 	for (i = 0; i < nitems(ld_pi_odd); i++) {
202 		ATF_CHECK(fabsl(sinl(ld_pi_odd[i])) < LDBL_EPSILON);
203 		ATF_CHECK(cosl(ld_pi_odd[i]) == -1.0);
204 		ATF_CHECK(fabsl(tanl(ld_pi_odd[i])) < LDBL_EPSILON);
205 
206 		ATF_CHECK(fabsl(sinl(-ld_pi_odd[i])) < LDBL_EPSILON);
207 		ATF_CHECK(cosl(-ld_pi_odd[i]) == -1.0);
208 		ATF_CHECK(fabsl(tanl(-ld_pi_odd[i])) < LDBL_EPSILON);
209 
210 		ATF_CHECK(fabsl(sinl(ld_pi_odd[i] * 2)) < LDBL_EPSILON);
211 		ATF_CHECK(cosl(ld_pi_odd[i] * 2) == 1.0);
212 		ATF_CHECK(fabsl(tanl(ld_pi_odd[i] * 2)) < LDBL_EPSILON);
213 
214 		ATF_CHECK(fabsl(sinl(-ld_pi_odd[i] * 2)) < LDBL_EPSILON);
215 		ATF_CHECK(cosl(-ld_pi_odd[i] * 2) == 1.0);
216 		ATF_CHECK(fabsl(tanl(-ld_pi_odd[i] * 2)) < LDBL_EPSILON);
217 	}
218 #endif
219 }
220 
221 ATF_TC(accuracy);
222 ATF_TC_HEAD(accuracy, tc)
223 {
224 
225 	atf_tc_set_md_var(tc, "descr",
226 	    "tests the accuracy of these functions over the primary range");
227 }
228 ATF_TC_BODY(accuracy, tc)
229 {
230 
231 	/* For small args, sin(x) = tan(x) = x, and cos(x) = 1. */
232 	testall(sin, 0xd.50ee515fe4aea16p-114L, 0xd.50ee515fe4aea16p-114L,
233 	     ALL_STD_EXCEPT, FE_INEXACT);
234 	testall(tan, 0xd.50ee515fe4aea16p-114L, 0xd.50ee515fe4aea16p-114L,
235 	     ALL_STD_EXCEPT, FE_INEXACT);
236 	testall(cos, 0xd.50ee515fe4aea16p-114L, 1.0,
237 		ALL_STD_EXCEPT, FE_INEXACT);
238 
239 	/*
240 	 * These tests should pass for f32, d64, and ld80 as long as
241 	 * the error is <= 0.75 ulp (round to nearest)
242 	 */
243 #if LDBL_MANT_DIG <= 64
244 #define	testacc	testall
245 #else
246 #define	testacc	testdf
247 #endif
248 	testacc(sin, 0.17255452780841205174L, 0.17169949801444412683L,
249 		ALL_STD_EXCEPT, FE_INEXACT);
250 	testacc(sin, -0.75431944555904520893L, -0.68479288156557286353L,
251 		ALL_STD_EXCEPT, FE_INEXACT);
252 	testacc(cos, 0.70556358769838947292L, 0.76124620693117771850L,
253 		ALL_STD_EXCEPT, FE_INEXACT);
254 	testacc(cos, -0.34061437849088045332L, 0.94254960031831729956L,
255 		ALL_STD_EXCEPT, FE_INEXACT);
256 	testacc(tan, -0.15862817413325692897L, -0.15997221861309522115L,
257 		ALL_STD_EXCEPT, FE_INEXACT);
258 	testacc(tan, 0.38374784931303813530L, 0.40376500259976759951L,
259 		ALL_STD_EXCEPT, FE_INEXACT);
260 
261 	/*
262 	 * XXX missing:
263 	 * - tests for ld128
264 	 * - tests for other rounding modes (probably won't pass for now)
265 	 * - tests for large numbers that get reduced to hi+lo with lo!=0
266 	 */
267 }
268 #endif
269 
270 ATF_TP_ADD_TCS(tp)
271 {
272 
273 	ATF_TP_ADD_TC(tp, special);
274 
275 #ifndef __i386__
276 	ATF_TP_ADD_TC(tp, accuracy);
277 	ATF_TP_ADD_TC(tp, reduction);
278 #endif
279 
280 	return (atf_no_error());
281 }
282