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 __FBSDID("$FreeBSD$");
38
39 #include <sys/param.h>
40
41 #include <assert.h>
42 #include <fenv.h>
43 #include <float.h>
44 #include <math.h>
45 #include <stdio.h>
46
47 #include <atf-c.h>
48
49 #include "test-utils.h"
50
51 #pragma STDC FENV_ACCESS ON
52
53 /*
54 * Test that a function returns the correct value and sets the
55 * exception flags correctly. The exceptmask specifies which
56 * exceptions we should check. We need to be lenient for several
57 * reasons, but mainly because on some architectures it's impossible
58 * to raise FE_OVERFLOW without raising FE_INEXACT.
59 *
60 * These are macros instead of functions so that assert provides more
61 * meaningful error messages.
62 *
63 * XXX The volatile here is to avoid gcc's bogus constant folding and work
64 * around the lack of support for the FENV_ACCESS pragma.
65 */
66 #define test(func, x, result, exceptmask, excepts) do { \
67 volatile long double _d = x; \
68 ATF_CHECK(feclearexcept(FE_ALL_EXCEPT) == 0); \
69 ATF_CHECK(fpequal((func)(_d), (result))); \
70 ATF_CHECK(((void)(func), fetestexcept(exceptmask) == (excepts))); \
71 } while (0)
72
73 #define testall(prefix, x, result, exceptmask, excepts) do { \
74 test(prefix, x, (double)result, exceptmask, excepts); \
75 test(prefix##f, x, (float)result, exceptmask, excepts); \
76 test(prefix##l, x, result, exceptmask, excepts); \
77 } while (0)
78
79 #define testdf(prefix, x, result, exceptmask, excepts) do { \
80 test(prefix, x, (double)result, exceptmask, excepts); \
81 test(prefix##f, x, (float)result, exceptmask, excepts); \
82 } while (0)
83
84 ATF_TC(special);
ATF_TC_HEAD(special,tc)85 ATF_TC_HEAD(special, tc)
86 {
87
88 atf_tc_set_md_var(tc, "descr",
89 "test special cases in sin(), cos(), and tan()");
90 }
ATF_TC_BODY(special,tc)91 ATF_TC_BODY(special, tc)
92 {
93
94 /* Values at 0 should be exact. */
95 testall(tan, 0.0, 0.0, ALL_STD_EXCEPT, 0);
96 testall(tan, -0.0, -0.0, ALL_STD_EXCEPT, 0);
97 testall(cos, 0.0, 1.0, ALL_STD_EXCEPT, 0);
98 testall(cos, -0.0, 1.0, ALL_STD_EXCEPT, 0);
99 testall(sin, 0.0, 0.0, ALL_STD_EXCEPT, 0);
100 testall(sin, -0.0, -0.0, ALL_STD_EXCEPT, 0);
101
102 /* func(+-Inf) == NaN */
103 testall(tan, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
104 testall(sin, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
105 testall(cos, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
106 testall(tan, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
107 testall(sin, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
108 testall(cos, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
109
110 /* func(NaN) == NaN */
111 testall(tan, NAN, NAN, ALL_STD_EXCEPT, 0);
112 testall(sin, NAN, NAN, ALL_STD_EXCEPT, 0);
113 testall(cos, NAN, NAN, ALL_STD_EXCEPT, 0);
114 }
115
116 #ifndef __i386__
117 ATF_TC(reduction);
ATF_TC_HEAD(reduction,tc)118 ATF_TC_HEAD(reduction, tc)
119 {
120
121 atf_tc_set_md_var(tc, "descr",
122 "tests to ensure argument reduction for large arguments is accurate");
123 }
ATF_TC_BODY(reduction,tc)124 ATF_TC_BODY(reduction, tc)
125 {
126 /* floats very close to odd multiples of pi */
127 static const float f_pi_odd[] = {
128 85563208.0f,
129 43998769152.0f,
130 9.2763667655669323e+25f,
131 1.5458357838905804e+29f,
132 };
133 /* doubles very close to odd multiples of pi */
134 static const double d_pi_odd[] = {
135 3.1415926535897931,
136 91.106186954104004,
137 642615.9188844458,
138 3397346.5699258847,
139 6134899525417045.0,
140 3.0213551960457761e+43,
141 1.2646209897993783e+295,
142 6.2083625380677099e+307,
143 };
144 /* long doubles very close to odd multiples of pi */
145 #if LDBL_MANT_DIG == 64
146 static const long double ld_pi_odd[] = {
147 1.1891886960373841596e+101L,
148 1.07999475322710967206e+2087L,
149 6.522151627890431836e+2147L,
150 8.9368974898260328229e+2484L,
151 9.2961044110572205863e+2555L,
152 4.90208421886578286e+3189L,
153 1.5275546401232615884e+3317L,
154 1.7227465626338900093e+3565L,
155 2.4160090594000745334e+3808L,
156 9.8477555741888350649e+4314L,
157 1.6061597222105160737e+4326L,
158 };
159 #endif
160
161 unsigned i;
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
220 ATF_TC(accuracy);
ATF_TC_HEAD(accuracy,tc)221 ATF_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 }
ATF_TC_BODY(accuracy,tc)227 ATF_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
ATF_TP_ADD_TCS(tp)269 ATF_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