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