1 /* $OpenBSD: trig_test.c,v 1.3 2021/12/13 18:04:28 deraadt Exp $ */
2 /*-
3 * Copyright (c) 2008 David Schultz <das@FreeBSD.org>
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above copyright
12 * notice, this list of conditions and the following disclaimer in the
13 * documentation and/or other materials provided with the distribution.
14 *
15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25 * SUCH DAMAGE.
26 */
27
28 #include "macros.h"
29
30 /*
31 * Tests for corner cases in trigonometric functions. Some accuracy tests
32 * are included as well, but these are very basic sanity checks, not
33 * intended to be comprehensive.
34 *
35 * The program for generating representable numbers near multiples of pi is
36 * available at http://www.cs.berkeley.edu/~wkahan/testpi/ .
37 */
38
39 #include <sys/types.h>
40
41 #include <fenv.h>
42 #include <float.h>
43 #include <math.h>
44 #include <stdio.h>
45
46 #include "test-utils.h"
47
48 #pragma STDC FENV_ACCESS ON
49
50 /*
51 * Test that a function returns the correct value and sets the
52 * exception flags correctly. The exceptmask specifies which
53 * exceptions we should check. We need to be lenient for several
54 * reasons, but mainly because on some architectures it's impossible
55 * to raise FE_OVERFLOW without raising FE_INEXACT.
56 *
57 * These are macros instead of functions so that assert provides more
58 * meaningful error messages.
59 *
60 * XXX The volatile here is to avoid gcc's bogus constant folding and work
61 * around the lack of support for the FENV_ACCESS pragma.
62 */
63 #define test(func, x, result, exceptmask, excepts) do { \
64 volatile long double _d = x; \
65 ATF_CHECK(feclearexcept(FE_ALL_EXCEPT) == 0); \
66 CHECK_FPEQUAL((func)(_d), (result)); \
67 CHECK_FP_EXCEPTIONS_MSG(excepts, exceptmask, "for %s(%s)", \
68 #func, #x); \
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 ATF_TC(special);
ATF_TC_HEAD(special,tc)83 ATF_TC_HEAD(special, tc)
84 {
85
86 atf_tc_set_md_var(tc, "descr",
87 "test special cases in sin(), cos(), and tan()");
88 }
ATF_TC_BODY(special,tc)89 ATF_TC_BODY(special, tc)
90 {
91
92 /* Values at 0 should be exact. */
93 testall(tan, 0.0, 0.0, ALL_STD_EXCEPT, 0);
94 testall(tan, -0.0, -0.0, ALL_STD_EXCEPT, 0);
95 testall(cos, 0.0, 1.0, ALL_STD_EXCEPT, 0);
96 testall(cos, -0.0, 1.0, ALL_STD_EXCEPT, 0);
97 testall(sin, 0.0, 0.0, ALL_STD_EXCEPT, 0);
98 testall(sin, -0.0, -0.0, ALL_STD_EXCEPT, 0);
99
100 /* func(+-Inf) == NaN */
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 testall(tan, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
105 testall(sin, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
106 testall(cos, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID);
107
108 /* func(NaN) == NaN */
109 testall(tan, NAN, NAN, ALL_STD_EXCEPT, 0);
110 testall(sin, NAN, NAN, ALL_STD_EXCEPT, 0);
111 testall(cos, NAN, NAN, ALL_STD_EXCEPT, 0);
112 }
113
114 #ifndef __i386__
115 ATF_TC(reduction);
ATF_TC_HEAD(reduction,tc)116 ATF_TC_HEAD(reduction, tc)
117 {
118
119 atf_tc_set_md_var(tc, "descr",
120 "tests to ensure argument reduction for large arguments is accurate");
121 }
ATF_TC_BODY(reduction,tc)122 ATF_TC_BODY(reduction, tc)
123 {
124 /* floats very close to odd multiples of pi */
125 static const float f_pi_odd[] = {
126 85563208.0f,
127 43998769152.0f,
128 9.2763667655669323e+25f,
129 1.5458357838905804e+29f,
130 };
131 /* doubles very close to odd multiples of pi */
132 static const double d_pi_odd[] = {
133 3.1415926535897931,
134 91.106186954104004,
135 642615.9188844458,
136 3397346.5699258847,
137 6134899525417045.0,
138 3.0213551960457761e+43,
139 1.2646209897993783e+295,
140 6.2083625380677099e+307,
141 };
142 /* long doubles very close to odd multiples of pi */
143 #if LDBL_MANT_DIG == 64
144 static const long double ld_pi_odd[] = {
145 1.1891886960373841596e+101L,
146 1.07999475322710967206e+2087L,
147 6.522151627890431836e+2147L,
148 8.9368974898260328229e+2484L,
149 9.2961044110572205863e+2555L,
150 4.90208421886578286e+3189L,
151 1.5275546401232615884e+3317L,
152 1.7227465626338900093e+3565L,
153 2.4160090594000745334e+3808L,
154 9.8477555741888350649e+4314L,
155 1.6061597222105160737e+4326L,
156 };
157 #endif
158
159 unsigned i;
160
161 #if defined(__amd64__) && defined(__clang__) && __clang_major__ >= 7 && \
162 __clang_major__ < 10 && __FreeBSD_cc_version < 1300002
163 atf_tc_expect_fail("test fails with clang 7-9 - bug 234040");
164 #endif
165
166 for (i = 0; i < nitems(f_pi_odd); i++) {
167 ATF_CHECK(fabs(sinf(f_pi_odd[i])) < FLT_EPSILON);
168 ATF_CHECK(cosf(f_pi_odd[i]) == -1.0);
169 ATF_CHECK(fabs(tan(f_pi_odd[i])) < FLT_EPSILON);
170
171 ATF_CHECK(fabs(sinf(-f_pi_odd[i])) < FLT_EPSILON);
172 ATF_CHECK(cosf(-f_pi_odd[i]) == -1.0);
173 ATF_CHECK(fabs(tanf(-f_pi_odd[i])) < FLT_EPSILON);
174
175 ATF_CHECK(fabs(sinf(f_pi_odd[i] * 2)) < FLT_EPSILON);
176 ATF_CHECK(cosf(f_pi_odd[i] * 2) == 1.0);
177 ATF_CHECK(fabs(tanf(f_pi_odd[i] * 2)) < FLT_EPSILON);
178
179 ATF_CHECK(fabs(sinf(-f_pi_odd[i] * 2)) < FLT_EPSILON);
180 ATF_CHECK(cosf(-f_pi_odd[i] * 2) == 1.0);
181 ATF_CHECK(fabs(tanf(-f_pi_odd[i] * 2)) < FLT_EPSILON);
182 }
183
184 for (i = 0; i < nitems(d_pi_odd); i++) {
185 ATF_CHECK(fabs(sin(d_pi_odd[i])) < 2 * DBL_EPSILON);
186 ATF_CHECK(cos(d_pi_odd[i]) == -1.0);
187 ATF_CHECK(fabs(tan(d_pi_odd[i])) < 2 * DBL_EPSILON);
188
189 ATF_CHECK(fabs(sin(-d_pi_odd[i])) < 2 * DBL_EPSILON);
190 ATF_CHECK(cos(-d_pi_odd[i]) == -1.0);
191 ATF_CHECK(fabs(tan(-d_pi_odd[i])) < 2 * DBL_EPSILON);
192
193 ATF_CHECK(fabs(sin(d_pi_odd[i] * 2)) < 2 * DBL_EPSILON);
194 ATF_CHECK(cos(d_pi_odd[i] * 2) == 1.0);
195 ATF_CHECK(fabs(tan(d_pi_odd[i] * 2)) < 2 * DBL_EPSILON);
196
197 ATF_CHECK(fabs(sin(-d_pi_odd[i] * 2)) < 2 * DBL_EPSILON);
198 ATF_CHECK(cos(-d_pi_odd[i] * 2) == 1.0);
199 ATF_CHECK(fabs(tan(-d_pi_odd[i] * 2)) < 2 * DBL_EPSILON);
200 }
201
202 #if LDBL_MANT_DIG == 64 /* XXX: || LDBL_MANT_DIG == 113 */
203 for (i = 0; i < nitems(ld_pi_odd); i++) {
204 ATF_CHECK(fabsl(sinl(ld_pi_odd[i])) < LDBL_EPSILON);
205 ATF_CHECK(cosl(ld_pi_odd[i]) == -1.0);
206 ATF_CHECK(fabsl(tanl(ld_pi_odd[i])) < LDBL_EPSILON);
207
208 ATF_CHECK(fabsl(sinl(-ld_pi_odd[i])) < LDBL_EPSILON);
209 ATF_CHECK(cosl(-ld_pi_odd[i]) == -1.0);
210 ATF_CHECK(fabsl(tanl(-ld_pi_odd[i])) < LDBL_EPSILON);
211
212 ATF_CHECK(fabsl(sinl(ld_pi_odd[i] * 2)) < LDBL_EPSILON);
213 ATF_CHECK(cosl(ld_pi_odd[i] * 2) == 1.0);
214 ATF_CHECK(fabsl(tanl(ld_pi_odd[i] * 2)) < LDBL_EPSILON);
215
216 ATF_CHECK(fabsl(sinl(-ld_pi_odd[i] * 2)) < LDBL_EPSILON);
217 ATF_CHECK(cosl(-ld_pi_odd[i] * 2) == 1.0);
218 ATF_CHECK(fabsl(tanl(-ld_pi_odd[i] * 2)) < LDBL_EPSILON);
219 }
220 #endif
221 }
222
223 ATF_TC(accuracy);
ATF_TC_HEAD(accuracy,tc)224 ATF_TC_HEAD(accuracy, tc)
225 {
226
227 atf_tc_set_md_var(tc, "descr",
228 "tests the accuracy of these functions over the primary range");
229 }
ATF_TC_BODY(accuracy,tc)230 ATF_TC_BODY(accuracy, tc)
231 {
232
233 /* For small args, sin(x) = tan(x) = x, and cos(x) = 1. */
234 testall(sin, 0xd.50ee515fe4aea16p-114L, 0xd.50ee515fe4aea16p-114L,
235 ALL_STD_EXCEPT, FE_INEXACT);
236 testall(tan, 0xd.50ee515fe4aea16p-114L, 0xd.50ee515fe4aea16p-114L,
237 ALL_STD_EXCEPT, FE_INEXACT);
238 testall(cos, 0xd.50ee515fe4aea16p-114L, 1.0,
239 ALL_STD_EXCEPT, FE_INEXACT);
240
241 /*
242 * These tests should pass for f32, d64, and ld80 as long as
243 * the error is <= 0.75 ulp (round to nearest)
244 */
245 #if LDBL_MANT_DIG <= 64
246 #define testacc testall
247 #else
248 #define testacc testdf
249 #endif
250 testacc(sin, 0.17255452780841205174L, 0.17169949801444412683L,
251 ALL_STD_EXCEPT, FE_INEXACT);
252 testacc(sin, -0.75431944555904520893L, -0.68479288156557286353L,
253 ALL_STD_EXCEPT, FE_INEXACT);
254 testacc(cos, 0.70556358769838947292L, 0.76124620693117771850L,
255 ALL_STD_EXCEPT, FE_INEXACT);
256 testacc(cos, -0.34061437849088045332L, 0.94254960031831729956L,
257 ALL_STD_EXCEPT, FE_INEXACT);
258 testacc(tan, -0.15862817413325692897L, -0.15997221861309522115L,
259 ALL_STD_EXCEPT, FE_INEXACT);
260 testacc(tan, 0.38374784931303813530L, 0.40376500259976759951L,
261 ALL_STD_EXCEPT, FE_INEXACT);
262
263 /*
264 * XXX missing:
265 * - tests for ld128
266 * - tests for other rounding modes (probably won't pass for now)
267 * - tests for large numbers that get reduced to hi+lo with lo!=0
268 */
269 }
270 #endif
271
ATF_TP_ADD_TCS(tp)272 ATF_TP_ADD_TCS(tp)
273 {
274
275 ATF_TP_ADD_TC(tp, special);
276
277 #ifndef __i386__
278 ATF_TP_ADD_TC(tp, accuracy);
279 ATF_TP_ADD_TC(tp, reduction);
280 #endif
281
282 return (atf_no_error());
283 }
284