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 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 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 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 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 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