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