invtrig_test.c revision 294243
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 the inverse trigonometric functions. Some
29 * accuracy tests are included as well, but these are very basic
30 * sanity checks, not intended to be comprehensive.
31 */
32
33#include <sys/cdefs.h>
34__FBSDID("$FreeBSD: stable/10/lib/msun/tests/invtrig_test.c 294243 2016-01-18 03:55:40Z ngie $");
35
36#include <assert.h>
37#include <fenv.h>
38#include <float.h>
39#include <math.h>
40#include <stdio.h>
41
42#include "test-utils.h"
43
44#pragma STDC FENV_ACCESS ON
45
46/*
47 * Test that a function returns the correct value and sets the
48 * exception flags correctly. A tolerance specifying the maximum
49 * relative error allowed may be specified. For the 'testall'
50 * functions, the tolerance is specified in ulps.
51 *
52 * These are macros instead of functions so that assert provides more
53 * meaningful error messages.
54 */
55#define	test_tol(func, x, result, tol, excepts) do {			\
56	volatile long double _in = (x), _out = (result);		\
57	assert(feclearexcept(FE_ALL_EXCEPT) == 0);			\
58	assert(fpequal_tol(func(_in), _out, (tol), CS_BOTH));		\
59	assert(((void)func, fetestexcept(ALL_STD_EXCEPT) == (excepts))); \
60} while (0)
61#define test(func, x, result, excepts)					\
62	test_tol(func, (x), (result), 0, (excepts))
63
64#define	_testall_tol(prefix, x, result, tol, excepts) do {		\
65	test_tol(prefix, (double)(x), (double)(result),			\
66		 (tol) * ldexp(1.0, 1 - DBL_MANT_DIG), (excepts));	\
67	test_tol(prefix##f, (float)(x), (float)(result),		\
68		 (tol) * ldexpf(1.0, 1 - FLT_MANT_DIG), (excepts));	\
69} while (0)
70
71#if LDBL_PREC == 53
72#define	testall_tol	_testall_tol
73#else
74#define	testall_tol(prefix, x, result, tol, excepts) do {		\
75	_testall_tol(prefix, x, result, tol, excepts);			\
76	test_tol(prefix##l, (x), (result),				\
77		 (tol) * ldexpl(1.0, 1 - LDBL_MANT_DIG), (excepts));	\
78} while (0)
79#endif
80
81#define testall(prefix, x, result, excepts)				\
82	testall_tol(prefix, (x), (result), 0, (excepts))
83
84#define	test2_tol(func, y, x, result, tol, excepts) do {		\
85	volatile long double _iny = (y), _inx = (x), _out = (result);	\
86	assert(feclearexcept(FE_ALL_EXCEPT) == 0);			\
87	assert(fpequal_tol(func(_iny, _inx), _out, (tol), CS_BOTH));	\
88	assert(((void)func, fetestexcept(ALL_STD_EXCEPT) == (excepts))); \
89} while (0)
90#define test2(func, y, x, result, excepts)				\
91	test2_tol(func, (y), (x), (result), 0, (excepts))
92
93#define	_testall2_tol(prefix, y, x, result, tol, excepts) do {		\
94	test2_tol(prefix, (double)(y), (double)(x), (double)(result),	\
95		  (tol) * ldexp(1.0, 1 - DBL_MANT_DIG), (excepts));	\
96	test2_tol(prefix##f, (float)(y), (float)(x), (float)(result),	\
97		  (tol) * ldexpf(1.0, 1 - FLT_MANT_DIG), (excepts));	\
98} while (0)
99
100#if LDBL_PREC == 53
101#define	testall2_tol	_testall2_tol
102#else
103#define	testall2_tol(prefix, y, x, result, tol, excepts) do {		\
104	_testall2_tol(prefix, y, x, result, tol, excepts);		\
105	test2_tol(prefix##l, (y), (x), (result),			\
106		  (tol) * ldexpl(1.0, 1 - LDBL_MANT_DIG), (excepts));	\
107} while (0)
108#endif
109
110#define testall2(prefix, y, x, result, excepts)				\
111	testall2_tol(prefix, (y), (x), (result), 0, (excepts))
112
113long double
114pi =   3.14159265358979323846264338327950280e+00L,
115pio3 = 1.04719755119659774615421446109316766e+00L,
116c3pi = 9.42477796076937971538793014983850839e+00L,
117c5pi = 1.57079632679489661923132169163975140e+01L,
118c7pi = 2.19911485751285526692385036829565196e+01L,
119c5pio3 = 5.23598775598298873077107230546583851e+00L,
120sqrt2m1 = 4.14213562373095048801688724209698081e-01L;
121
122
123/*
124 * Test special case inputs in asin(), acos() and atan(): signed
125 * zeroes, infinities, and NaNs.
126 */
127static void
128test_special(void)
129{
130
131	testall(asin, 0.0, 0.0, 0);
132	testall(acos, 0.0, pi / 2, FE_INEXACT);
133	testall(atan, 0.0, 0.0, 0);
134	testall(asin, -0.0, -0.0, 0);
135	testall(acos, -0.0, pi / 2, FE_INEXACT);
136	testall(atan, -0.0, -0.0, 0);
137
138	testall(asin, INFINITY, NAN, FE_INVALID);
139	testall(acos, INFINITY, NAN, FE_INVALID);
140	testall(atan, INFINITY, pi / 2, FE_INEXACT);
141	testall(asin, -INFINITY, NAN, FE_INVALID);
142	testall(acos, -INFINITY, NAN, FE_INVALID);
143	testall(atan, -INFINITY, -pi / 2, FE_INEXACT);
144
145	testall(asin, NAN, NAN, 0);
146	testall(acos, NAN, NAN, 0);
147	testall(atan, NAN, NAN, 0);
148}
149
150/*
151 * Test special case inputs in atan2(), where the exact value of y/x is
152 * zero or non-finite.
153 */
154static void
155test_special_atan2(void)
156{
157	long double z;
158	int e;
159
160	testall2(atan2, 0.0, -0.0, pi, FE_INEXACT);
161	testall2(atan2, -0.0, -0.0, -pi, FE_INEXACT);
162	testall2(atan2, 0.0, 0.0, 0.0, 0);
163	testall2(atan2, -0.0, 0.0, -0.0, 0);
164
165	testall2(atan2, INFINITY, -INFINITY, c3pi / 4, FE_INEXACT);
166	testall2(atan2, -INFINITY, -INFINITY, -c3pi / 4, FE_INEXACT);
167	testall2(atan2, INFINITY, INFINITY, pi / 4, FE_INEXACT);
168	testall2(atan2, -INFINITY, INFINITY, -pi / 4, FE_INEXACT);
169
170	/* Tests with one input in the range (0, Inf]. */
171	z = 1.23456789L;
172	for (e = FLT_MIN_EXP - FLT_MANT_DIG; e <= FLT_MAX_EXP; e++) {
173		test2(atan2f, 0.0, ldexpf(z, e), 0.0, 0);
174		test2(atan2f, -0.0, ldexpf(z, e), -0.0, 0);
175		test2(atan2f, 0.0, ldexpf(-z, e), (float)pi, FE_INEXACT);
176		test2(atan2f, -0.0, ldexpf(-z, e), (float)-pi, FE_INEXACT);
177		test2(atan2f, ldexpf(z, e), 0.0, (float)pi / 2, FE_INEXACT);
178		test2(atan2f, ldexpf(z, e), -0.0, (float)pi / 2, FE_INEXACT);
179		test2(atan2f, ldexpf(-z, e), 0.0, (float)-pi / 2, FE_INEXACT);
180		test2(atan2f, ldexpf(-z, e), -0.0, (float)-pi / 2, FE_INEXACT);
181	}
182	for (e = DBL_MIN_EXP - DBL_MANT_DIG; e <= DBL_MAX_EXP; e++) {
183		test2(atan2, 0.0, ldexp(z, e), 0.0, 0);
184		test2(atan2, -0.0, ldexp(z, e), -0.0, 0);
185		test2(atan2, 0.0, ldexp(-z, e), (double)pi, FE_INEXACT);
186		test2(atan2, -0.0, ldexp(-z, e), (double)-pi, FE_INEXACT);
187		test2(atan2, ldexp(z, e), 0.0, (double)pi / 2, FE_INEXACT);
188		test2(atan2, ldexp(z, e), -0.0, (double)pi / 2, FE_INEXACT);
189		test2(atan2, ldexp(-z, e), 0.0, (double)-pi / 2, FE_INEXACT);
190		test2(atan2, ldexp(-z, e), -0.0, (double)-pi / 2, FE_INEXACT);
191	}
192	for (e = LDBL_MIN_EXP - LDBL_MANT_DIG; e <= LDBL_MAX_EXP; e++) {
193		test2(atan2l, 0.0, ldexpl(z, e), 0.0, 0);
194		test2(atan2l, -0.0, ldexpl(z, e), -0.0, 0);
195		test2(atan2l, 0.0, ldexpl(-z, e), pi, FE_INEXACT);
196		test2(atan2l, -0.0, ldexpl(-z, e), -pi, FE_INEXACT);
197		test2(atan2l, ldexpl(z, e), 0.0, pi / 2, FE_INEXACT);
198		test2(atan2l, ldexpl(z, e), -0.0, pi / 2, FE_INEXACT);
199		test2(atan2l, ldexpl(-z, e), 0.0, -pi / 2, FE_INEXACT);
200		test2(atan2l, ldexpl(-z, e), -0.0, -pi / 2, FE_INEXACT);
201	}
202
203	/* Tests with one input in the range (0, Inf). */
204	for (e = FLT_MIN_EXP - FLT_MANT_DIG; e <= FLT_MAX_EXP - 1; e++) {
205		test2(atan2f, ldexpf(z, e), INFINITY, 0.0, 0);
206		test2(atan2f, ldexpf(-z,e), INFINITY, -0.0, 0);
207		test2(atan2f, ldexpf(z, e), -INFINITY, (float)pi, FE_INEXACT);
208		test2(atan2f, ldexpf(-z,e), -INFINITY, (float)-pi, FE_INEXACT);
209		test2(atan2f, INFINITY, ldexpf(z,e), (float)pi/2, FE_INEXACT);
210		test2(atan2f, INFINITY, ldexpf(-z,e), (float)pi/2, FE_INEXACT);
211		test2(atan2f, -INFINITY, ldexpf(z,e), (float)-pi/2,FE_INEXACT);
212		test2(atan2f, -INFINITY, ldexpf(-z,e),(float)-pi/2,FE_INEXACT);
213	}
214	for (e = DBL_MIN_EXP - DBL_MANT_DIG; e <= DBL_MAX_EXP - 1; e++) {
215		test2(atan2, ldexp(z, e), INFINITY, 0.0, 0);
216		test2(atan2, ldexp(-z,e), INFINITY, -0.0, 0);
217		test2(atan2, ldexp(z, e), -INFINITY, (double)pi, FE_INEXACT);
218		test2(atan2, ldexp(-z,e), -INFINITY, (double)-pi, FE_INEXACT);
219		test2(atan2, INFINITY, ldexp(z,e), (double)pi/2, FE_INEXACT);
220		test2(atan2, INFINITY, ldexp(-z,e), (double)pi/2, FE_INEXACT);
221		test2(atan2, -INFINITY, ldexp(z,e), (double)-pi/2,FE_INEXACT);
222		test2(atan2, -INFINITY, ldexp(-z,e),(double)-pi/2,FE_INEXACT);
223	}
224	for (e = LDBL_MIN_EXP - LDBL_MANT_DIG; e <= LDBL_MAX_EXP - 1; e++) {
225		test2(atan2l, ldexpl(z, e), INFINITY, 0.0, 0);
226		test2(atan2l, ldexpl(-z,e), INFINITY, -0.0, 0);
227		test2(atan2l, ldexpl(z, e), -INFINITY, pi, FE_INEXACT);
228		test2(atan2l, ldexpl(-z,e), -INFINITY, -pi, FE_INEXACT);
229		test2(atan2l, INFINITY, ldexpl(z, e), pi / 2, FE_INEXACT);
230		test2(atan2l, INFINITY, ldexpl(-z, e), pi / 2, FE_INEXACT);
231		test2(atan2l, -INFINITY, ldexpl(z, e), -pi / 2, FE_INEXACT);
232		test2(atan2l, -INFINITY, ldexpl(-z, e), -pi / 2, FE_INEXACT);
233	}
234}
235
236/*
237 * Test various inputs to asin(), acos() and atan() and verify that the
238 * results are accurate to within 1 ulp.
239 */
240static void
241test_accuracy(void)
242{
243
244	/* We expect correctly rounded results for these basic cases. */
245	testall(asin, 1.0, pi / 2, FE_INEXACT);
246	testall(acos, 1.0, 0, 0);
247	testall(atan, 1.0, pi / 4, FE_INEXACT);
248	testall(asin, -1.0, -pi / 2, FE_INEXACT);
249	testall(acos, -1.0, pi, FE_INEXACT);
250	testall(atan, -1.0, -pi / 4, FE_INEXACT);
251
252	/*
253	 * Here we expect answers to be within 1 ulp, although inexactness
254	 * in the input, combined with double rounding, could cause larger
255	 * errors.
256	 */
257
258	testall_tol(asin, sqrtl(2) / 2, pi / 4, 1, FE_INEXACT);
259	testall_tol(acos, sqrtl(2) / 2, pi / 4, 1, FE_INEXACT);
260	testall_tol(asin, -sqrtl(2) / 2, -pi / 4, 1, FE_INEXACT);
261	testall_tol(acos, -sqrtl(2) / 2, c3pi / 4, 1, FE_INEXACT);
262
263	testall_tol(asin, sqrtl(3) / 2, pio3, 1, FE_INEXACT);
264	testall_tol(acos, sqrtl(3) / 2, pio3 / 2, 1, FE_INEXACT);
265	testall_tol(atan, sqrtl(3), pio3, 1, FE_INEXACT);
266	testall_tol(asin, -sqrtl(3) / 2, -pio3, 1, FE_INEXACT);
267	testall_tol(acos, -sqrtl(3) / 2, c5pio3 / 2, 1, FE_INEXACT);
268	testall_tol(atan, -sqrtl(3), -pio3, 1, FE_INEXACT);
269
270	testall_tol(atan, sqrt2m1, pi / 8, 1, FE_INEXACT);
271	testall_tol(atan, -sqrt2m1, -pi / 8, 1, FE_INEXACT);
272}
273
274/*
275 * Test inputs to atan2() where x is a power of 2. These are easy cases
276 * because y/x is exact.
277 */
278static void
279test_p2x_atan2(void)
280{
281
282	testall2(atan2, 1.0, 1.0, pi / 4, FE_INEXACT);
283	testall2(atan2, 1.0, -1.0, c3pi / 4, FE_INEXACT);
284	testall2(atan2, -1.0, 1.0, -pi / 4, FE_INEXACT);
285	testall2(atan2, -1.0, -1.0, -c3pi / 4, FE_INEXACT);
286
287	testall2_tol(atan2, sqrt2m1 * 2, 2.0, pi / 8, 1, FE_INEXACT);
288	testall2_tol(atan2, sqrt2m1 * 2, -2.0, c7pi / 8, 1, FE_INEXACT);
289	testall2_tol(atan2, -sqrt2m1 * 2, 2.0, -pi / 8, 1, FE_INEXACT);
290	testall2_tol(atan2, -sqrt2m1 * 2, -2.0, -c7pi / 8, 1, FE_INEXACT);
291
292	testall2_tol(atan2, sqrtl(3) * 0.5, 0.5, pio3, 1, FE_INEXACT);
293	testall2_tol(atan2, sqrtl(3) * 0.5, -0.5, pio3 * 2, 1, FE_INEXACT);
294	testall2_tol(atan2, -sqrtl(3) * 0.5, 0.5, -pio3, 1, FE_INEXACT);
295	testall2_tol(atan2, -sqrtl(3) * 0.5, -0.5, -pio3 * 2, 1, FE_INEXACT);
296}
297
298/*
299 * Test inputs very close to 0.
300 */
301static void
302test_tiny(void)
303{
304	float tiny = 0x1.23456p-120f;
305
306	testall(asin, tiny, tiny, FE_INEXACT);
307	testall(acos, tiny, pi / 2, FE_INEXACT);
308	testall(atan, tiny, tiny, FE_INEXACT);
309
310	testall(asin, -tiny, -tiny, FE_INEXACT);
311	testall(acos, -tiny, pi / 2, FE_INEXACT);
312	testall(atan, -tiny, -tiny, FE_INEXACT);
313
314	/* Test inputs to atan2() that would cause y/x to underflow. */
315	test2(atan2f, 0x1.0p-100, 0x1.0p100, 0.0, FE_INEXACT | FE_UNDERFLOW);
316	test2(atan2, 0x1.0p-1000, 0x1.0p1000, 0.0, FE_INEXACT | FE_UNDERFLOW);
317	test2(atan2l, ldexpl(1.0, 100 - LDBL_MAX_EXP),
318	      ldexpl(1.0, LDBL_MAX_EXP - 100), 0.0, FE_INEXACT | FE_UNDERFLOW);
319	test2(atan2f, -0x1.0p-100, 0x1.0p100, -0.0, FE_INEXACT | FE_UNDERFLOW);
320	test2(atan2, -0x1.0p-1000, 0x1.0p1000, -0.0, FE_INEXACT | FE_UNDERFLOW);
321	test2(atan2l, -ldexpl(1.0, 100 - LDBL_MAX_EXP),
322	      ldexpl(1.0, LDBL_MAX_EXP - 100), -0.0, FE_INEXACT | FE_UNDERFLOW);
323	test2(atan2f, 0x1.0p-100, -0x1.0p100, (float)pi, FE_INEXACT);
324	test2(atan2, 0x1.0p-1000, -0x1.0p1000, (double)pi, FE_INEXACT);
325	test2(atan2l, ldexpl(1.0, 100 - LDBL_MAX_EXP),
326	      -ldexpl(1.0, LDBL_MAX_EXP - 100), pi, FE_INEXACT);
327	test2(atan2f, -0x1.0p-100, -0x1.0p100, (float)-pi, FE_INEXACT);
328	test2(atan2, -0x1.0p-1000, -0x1.0p1000, (double)-pi, FE_INEXACT);
329	test2(atan2l, -ldexpl(1.0, 100 - LDBL_MAX_EXP),
330	      -ldexpl(1.0, LDBL_MAX_EXP - 100), -pi, FE_INEXACT);
331}
332
333/*
334 * Test very large inputs to atan().
335 */
336static void
337test_atan_huge(void)
338{
339	float huge = 0x1.23456p120;
340
341	testall(atan, huge, pi / 2, FE_INEXACT);
342	testall(atan, -huge, -pi / 2, FE_INEXACT);
343
344	/* Test inputs to atan2() that would cause y/x to overflow. */
345	test2(atan2f, 0x1.0p100, 0x1.0p-100, (float)pi / 2, FE_INEXACT);
346	test2(atan2, 0x1.0p1000, 0x1.0p-1000, (double)pi / 2, FE_INEXACT);
347	test2(atan2l, ldexpl(1.0, LDBL_MAX_EXP - 100),
348	      ldexpl(1.0, 100 - LDBL_MAX_EXP), pi / 2, FE_INEXACT);
349	test2(atan2f, -0x1.0p100, 0x1.0p-100, (float)-pi / 2, FE_INEXACT);
350	test2(atan2, -0x1.0p1000, 0x1.0p-1000, (double)-pi / 2, FE_INEXACT);
351	test2(atan2l, -ldexpl(1.0, LDBL_MAX_EXP - 100),
352	      ldexpl(1.0, 100 - LDBL_MAX_EXP), -pi / 2, FE_INEXACT);
353
354	test2(atan2f, 0x1.0p100, -0x1.0p-100, (float)pi / 2, FE_INEXACT);
355	test2(atan2, 0x1.0p1000, -0x1.0p-1000, (double)pi / 2, FE_INEXACT);
356	test2(atan2l, ldexpl(1.0, LDBL_MAX_EXP - 100),
357	      -ldexpl(1.0, 100 - LDBL_MAX_EXP), pi / 2, FE_INEXACT);
358	test2(atan2f, -0x1.0p100, -0x1.0p-100, (float)-pi / 2, FE_INEXACT);
359	test2(atan2, -0x1.0p1000, -0x1.0p-1000, (double)-pi / 2, FE_INEXACT);
360	test2(atan2l, -ldexpl(1.0, LDBL_MAX_EXP - 100),
361	      -ldexpl(1.0, 100 - LDBL_MAX_EXP), -pi / 2, FE_INEXACT);
362}
363
364/*
365 * Test that sin(asin(x)) == x, and similarly for acos() and atan().
366 * You need to have a working sinl(), cosl(), and tanl() for these
367 * tests to pass.
368 */
369static long double
370sinasinf(float x)
371{
372
373	return (sinl(asinf(x)));
374}
375
376static long double
377sinasin(double x)
378{
379
380	return (sinl(asin(x)));
381}
382
383static long double
384sinasinl(long double x)
385{
386
387	return (sinl(asinl(x)));
388}
389
390static long double
391cosacosf(float x)
392{
393
394	return (cosl(acosf(x)));
395}
396
397static long double
398cosacos(double x)
399{
400
401	return (cosl(acos(x)));
402}
403
404static long double
405cosacosl(long double x)
406{
407
408	return (cosl(acosl(x)));
409}
410
411static long double
412tanatanf(float x)
413{
414
415	return (tanl(atanf(x)));
416}
417
418static long double
419tanatan(double x)
420{
421
422	return (tanl(atan(x)));
423}
424
425static long double
426tanatanl(long double x)
427{
428
429	return (tanl(atanl(x)));
430}
431
432static void
433test_inverse(void)
434{
435	float i;
436
437	for (i = -1; i <= 1; i += 0x1.0p-12f) {
438		testall_tol(sinasin, i, i, 2, i == 0 ? 0 : FE_INEXACT);
439		/* The relative error for cosacos is very large near x=0. */
440		if (fabsf(i) > 0x1.0p-4f)
441			testall_tol(cosacos, i, i, 16, i == 1 ? 0 : FE_INEXACT);
442		testall_tol(tanatan, i, i, 2, i == 0 ? 0 : FE_INEXACT);
443	}
444}
445
446int
447main(int argc, char *argv[])
448{
449
450#if defined(__i386__)
451	printf("1..0 # SKIP fails all assertions on i386\n");
452	return (0);
453#endif
454
455	printf("1..7\n");
456
457	test_special();
458	printf("ok 1 - special\n");
459
460	test_special_atan2();
461	printf("ok 2 - atan2 special\n");
462
463	test_accuracy();
464	printf("ok 3 - accuracy\n");
465
466	test_p2x_atan2();
467	printf("ok 4 - atan2 p2x\n");
468
469	test_tiny();
470	printf("ok 5 - tiny inputs\n");
471
472	test_atan_huge();
473	printf("ok 6 - atan huge inputs\n");
474
475	test_inverse();
476	printf("ok 7 - inverse\n");
477
478	return (0);
479}
480