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 82ATF_TC(special); 83ATF_TC_HEAD(special, tc) 84{ 85 86 atf_tc_set_md_var(tc, "descr", 87 "test special cases in sin(), cos(), and tan()"); 88} 89ATF_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__ 115ATF_TC(reduction); 116ATF_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} 122ATF_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 223ATF_TC(accuracy); 224ATF_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} 230ATF_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 272ATF_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