trig_test.c revision 176748
1176379Sdas/*- 2176379Sdas * Copyright (c) 2008 David Schultz <das@FreeBSD.org> 3176379Sdas * All rights reserved. 4176379Sdas * 5176379Sdas * Redistribution and use in source and binary forms, with or without 6176379Sdas * modification, are permitted provided that the following conditions 7176379Sdas * are met: 8176379Sdas * 1. Redistributions of source code must retain the above copyright 9176379Sdas * notice, this list of conditions and the following disclaimer. 10176379Sdas * 2. Redistributions in binary form must reproduce the above copyright 11176379Sdas * notice, this list of conditions and the following disclaimer in the 12176379Sdas * documentation and/or other materials provided with the distribution. 13176379Sdas * 14176379Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15176379Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16176379Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17176379Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18176379Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19176379Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20176379Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21176379Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22176379Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23176379Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24176379Sdas * SUCH DAMAGE. 25176379Sdas */ 26176379Sdas 27176379Sdas/* 28176379Sdas * Tests for corner cases in trigonometric functions. Some accuracy tests 29176379Sdas * are included as well, but these are very basic sanity checks, not 30176379Sdas * intended to be comprehensive. 31176379Sdas * 32176379Sdas * The program for generating representable numbers near multiples of pi is 33176379Sdas * available at http://www.cs.berkeley.edu/~wkahan/testpi/ . 34176379Sdas */ 35176379Sdas 36176379Sdas#include <sys/cdefs.h> 37176379Sdas__FBSDID("$FreeBSD: head/tools/regression/lib/msun/test-trig.c 176748 2008-03-02 20:49:24Z das $"); 38176379Sdas 39176379Sdas#include <assert.h> 40176379Sdas#include <fenv.h> 41176379Sdas#include <float.h> 42176379Sdas#include <math.h> 43176379Sdas#include <stdio.h> 44176379Sdas 45176379Sdas#define ALL_STD_EXCEPT (FE_DIVBYZERO | FE_INEXACT | FE_INVALID | \ 46176379Sdas FE_OVERFLOW | FE_UNDERFLOW) 47176379Sdas 48176379Sdas#define LEN(a) (sizeof(a) / sizeof((a)[0])) 49176379Sdas 50176379Sdas#pragma STDC FENV_ACCESS ON 51176379Sdas 52176379Sdas/* 53176379Sdas * Test that a function returns the correct value and sets the 54176379Sdas * exception flags correctly. The exceptmask specifies which 55176379Sdas * exceptions we should check. We need to be lenient for several 56176379Sdas * reasons, but mainly because on some architectures it's impossible 57176379Sdas * to raise FE_OVERFLOW without raising FE_INEXACT. 58176379Sdas * 59176379Sdas * These are macros instead of functions so that assert provides more 60176379Sdas * meaningful error messages. 61176379Sdas * 62176379Sdas * XXX The volatile here is to avoid gcc's bogus constant folding and work 63176379Sdas * around the lack of support for the FENV_ACCESS pragma. 64176379Sdas */ 65176379Sdas#define test(func, x, result, exceptmask, excepts) do { \ 66176379Sdas volatile long double _d = x; \ 67176379Sdas assert(feclearexcept(FE_ALL_EXCEPT) == 0); \ 68176748Sdas assert(fpequal((func)(_d), (result))); \ 69176379Sdas assert(((func), fetestexcept(exceptmask) == (excepts))); \ 70176379Sdas} while (0) 71176379Sdas 72176379Sdas#define testall(prefix, x, result, exceptmask, excepts) do { \ 73176748Sdas test(prefix, x, (double)result, exceptmask, excepts); \ 74176379Sdas test(prefix##f, x, (float)result, exceptmask, excepts); \ 75176379Sdas test(prefix##l, x, result, exceptmask, excepts); \ 76176379Sdas} while (0) 77176379Sdas 78176748Sdas#define testdf(prefix, x, result, exceptmask, excepts) do { \ 79176748Sdas test(prefix, x, (double)result, exceptmask, excepts); \ 80176748Sdas test(prefix##f, x, (float)result, exceptmask, excepts); \ 81176748Sdas} while (0) 82176748Sdas 83176748Sdas 84176748Sdas 85176379Sdas/* 86176379Sdas * Determine whether x and y are equal, with two special rules: 87176379Sdas * +0.0 != -0.0 88176379Sdas * NaN == NaN 89176379Sdas */ 90176379Sdasint 91176379Sdasfpequal(long double x, long double y) 92176379Sdas{ 93176379Sdas return ((x == y && signbit(x) == signbit(y)) || isnan(x) && isnan(y)); 94176379Sdas} 95176379Sdas 96176379Sdas/* 97176379Sdas * Test special cases in sin(), cos(), and tan(). 98176379Sdas */ 99176379Sdasstatic void 100176379Sdasrun_special_tests(void) 101176379Sdas{ 102176379Sdas 103176379Sdas /* Values at 0 should be exact. */ 104176379Sdas testall(tan, 0.0, 0.0, ALL_STD_EXCEPT, 0); 105176379Sdas testall(tan, -0.0, -0.0, ALL_STD_EXCEPT, 0); 106176379Sdas testall(cos, 0.0, 1.0, ALL_STD_EXCEPT, 0); 107176379Sdas testall(cos, -0.0, 1.0, ALL_STD_EXCEPT, 0); 108176379Sdas testall(sin, 0.0, 0.0, ALL_STD_EXCEPT, 0); 109176379Sdas testall(sin, -0.0, -0.0, ALL_STD_EXCEPT, 0); 110176379Sdas 111176379Sdas /* func(+-Inf) == NaN */ 112176379Sdas testall(tan, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID); 113176379Sdas testall(sin, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID); 114176379Sdas testall(cos, INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID); 115176379Sdas testall(tan, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID); 116176379Sdas testall(sin, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID); 117176379Sdas testall(cos, -INFINITY, NAN, ALL_STD_EXCEPT, FE_INVALID); 118176379Sdas 119176379Sdas /* func(NaN) == NaN */ 120176379Sdas testall(tan, NAN, NAN, ALL_STD_EXCEPT, 0); 121176379Sdas testall(sin, NAN, NAN, ALL_STD_EXCEPT, 0); 122176379Sdas testall(cos, NAN, NAN, ALL_STD_EXCEPT, 0); 123176379Sdas} 124176379Sdas 125176379Sdas/* 126176379Sdas * Tests to ensure argument reduction for large arguments is accurate. 127176379Sdas */ 128176379Sdasstatic void 129176379Sdasrun_reduction_tests(void) 130176379Sdas{ 131176379Sdas /* floats very close to odd multiples of pi */ 132176379Sdas static const float f_pi_odd[] = { 133176379Sdas 85563208.0f, 134176379Sdas 43998769152.0f, 135176379Sdas 9.2763667655669323e+25f, 136176379Sdas 1.5458357838905804e+29f, 137176379Sdas }; 138176379Sdas /* doubles very close to odd multiples of pi */ 139176379Sdas static const double d_pi_odd[] = { 140176379Sdas 3.1415926535897931, 141176379Sdas 91.106186954104004, 142176379Sdas 642615.9188844458, 143176379Sdas 3397346.5699258847, 144176379Sdas 6134899525417045.0, 145176379Sdas 3.0213551960457761e+43, 146176379Sdas 1.2646209897993783e+295, 147176379Sdas 6.2083625380677099e+307, 148176379Sdas }; 149176379Sdas /* long doubles very close to odd multiples of pi */ 150176379Sdas#if LDBL_MANT_DIG == 64 151176379Sdas static const long double ld_pi_odd[] = { 152176379Sdas 1.1891886960373841596e+101L, 153176379Sdas 1.07999475322710967206e+2087L, 154176379Sdas 6.522151627890431836e+2147L, 155176379Sdas 8.9368974898260328229e+2484L, 156176379Sdas 9.2961044110572205863e+2555L, 157176379Sdas 4.90208421886578286e+3189L, 158176379Sdas 1.5275546401232615884e+3317L, 159176379Sdas 1.7227465626338900093e+3565L, 160176379Sdas 2.4160090594000745334e+3808L, 161176379Sdas 9.8477555741888350649e+4314L, 162176379Sdas 1.6061597222105160737e+4326L, 163176379Sdas }; 164176379Sdas#elif LDBL_MANT_DIG == 113 165176379Sdas static const long double ld_pi_odd[] = { 166176379Sdas /* XXX */ 167176379Sdas }; 168176379Sdas#endif 169176379Sdas 170176379Sdas int i; 171176379Sdas 172176379Sdas for (i = 0; i < LEN(f_pi_odd); i++) { 173176379Sdas assert(fabs(sinf(f_pi_odd[i])) < FLT_EPSILON); 174176379Sdas assert(cosf(f_pi_odd[i]) == -1.0); 175176379Sdas assert(fabs(tan(f_pi_odd[i])) < FLT_EPSILON); 176176379Sdas 177176379Sdas assert(fabs(sinf(-f_pi_odd[i])) < FLT_EPSILON); 178176379Sdas assert(cosf(-f_pi_odd[i]) == -1.0); 179176379Sdas assert(fabs(tanf(-f_pi_odd[i])) < FLT_EPSILON); 180176379Sdas 181176379Sdas assert(fabs(sinf(f_pi_odd[i] * 2)) < FLT_EPSILON); 182176379Sdas assert(cosf(f_pi_odd[i] * 2) == 1.0); 183176379Sdas assert(fabs(tanf(f_pi_odd[i] * 2)) < FLT_EPSILON); 184176379Sdas 185176379Sdas assert(fabs(sinf(-f_pi_odd[i] * 2)) < FLT_EPSILON); 186176379Sdas assert(cosf(-f_pi_odd[i] * 2) == 1.0); 187176379Sdas assert(fabs(tanf(-f_pi_odd[i] * 2)) < FLT_EPSILON); 188176379Sdas } 189176379Sdas 190176379Sdas for (i = 0; i < LEN(d_pi_odd); i++) { 191176379Sdas assert(fabs(sin(d_pi_odd[i])) < 2 * DBL_EPSILON); 192176379Sdas assert(cos(d_pi_odd[i]) == -1.0); 193176379Sdas assert(fabs(tan(d_pi_odd[i])) < 2 * DBL_EPSILON); 194176379Sdas 195176379Sdas assert(fabs(sin(-d_pi_odd[i])) < 2 * DBL_EPSILON); 196176379Sdas assert(cos(-d_pi_odd[i]) == -1.0); 197176379Sdas assert(fabs(tan(-d_pi_odd[i])) < 2 * DBL_EPSILON); 198176379Sdas 199176379Sdas assert(fabs(sin(d_pi_odd[i] * 2)) < 2 * DBL_EPSILON); 200176379Sdas assert(cos(d_pi_odd[i] * 2) == 1.0); 201176379Sdas assert(fabs(tan(d_pi_odd[i] * 2)) < 2 * DBL_EPSILON); 202176379Sdas 203176379Sdas assert(fabs(sin(-d_pi_odd[i] * 2)) < 2 * DBL_EPSILON); 204176379Sdas assert(cos(-d_pi_odd[i] * 2) == 1.0); 205176379Sdas assert(fabs(tan(-d_pi_odd[i] * 2)) < 2 * DBL_EPSILON); 206176379Sdas } 207176379Sdas 208176379Sdas#if LDBL_MANT_DIG > 53 209176379Sdas for (i = 0; i < LEN(ld_pi_odd); i++) { 210176379Sdas assert(fabsl(sinl(ld_pi_odd[i])) < LDBL_EPSILON); 211176379Sdas assert(cosl(ld_pi_odd[i]) == -1.0); 212176379Sdas assert(fabsl(tanl(ld_pi_odd[i])) < LDBL_EPSILON); 213176379Sdas 214176379Sdas assert(fabsl(sinl(-ld_pi_odd[i])) < LDBL_EPSILON); 215176379Sdas assert(cosl(-ld_pi_odd[i]) == -1.0); 216176379Sdas assert(fabsl(tanl(-ld_pi_odd[i])) < LDBL_EPSILON); 217176379Sdas 218176379Sdas assert(fabsl(sinl(ld_pi_odd[i] * 2)) < LDBL_EPSILON); 219176379Sdas assert(cosl(ld_pi_odd[i] * 2) == 1.0); 220176379Sdas assert(fabsl(tanl(ld_pi_odd[i] * 2)) < LDBL_EPSILON); 221176379Sdas 222176379Sdas assert(fabsl(sinl(-ld_pi_odd[i] * 2)) < LDBL_EPSILON); 223176379Sdas assert(cosl(-ld_pi_odd[i] * 2) == 1.0); 224176379Sdas assert(fabsl(tanl(-ld_pi_odd[i] * 2)) < LDBL_EPSILON); 225176379Sdas } 226176379Sdas#endif 227176379Sdas} 228176379Sdas 229176379Sdas/* 230176379Sdas * Tests the accuracy of these functions over the primary range. 231176379Sdas */ 232176379Sdasstatic void 233176379Sdasrun_accuracy_tests(void) 234176379Sdas{ 235176379Sdas 236176379Sdas /* For small args, sin(x) = tan(x) = x, and cos(x) = 1. */ 237176379Sdas testall(sin, 0xd.50ee515fe4aea16p-114L, 0xd.50ee515fe4aea16p-114L, 238176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 239176379Sdas testall(tan, 0xd.50ee515fe4aea16p-114L, 0xd.50ee515fe4aea16p-114L, 240176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 241176379Sdas testall(cos, 0xd.50ee515fe4aea16p-114L, 1.0, 242176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 243176379Sdas 244176379Sdas /* 245176379Sdas * These tests should pass for f32, d64, and ld80 as long as 246176379Sdas * the error is <= 0.75 ulp (round to nearest) 247176379Sdas */ 248176748Sdas#if LDBL_MANT_DIG <= 64 249176748Sdas#define testacc testall 250176748Sdas#else 251176748Sdas#define testacc testdf 252176748Sdas#endif 253176748Sdas testacc(sin, 0.17255452780841205174L, 0.17169949801444412683L, 254176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 255176748Sdas testacc(sin, -0.75431944555904520893L, -0.68479288156557286353L, 256176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 257176748Sdas testacc(cos, 0.70556358769838947292L, 0.76124620693117771850L, 258176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 259176748Sdas testacc(cos, -0.34061437849088045332L, 0.94254960031831729956L, 260176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 261176748Sdas testacc(tan, -0.15862817413325692897L, -0.15997221861309522115L, 262176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 263176748Sdas testacc(tan, 0.38374784931303813530L, 0.40376500259976759951L, 264176379Sdas ALL_STD_EXCEPT, FE_INEXACT); 265176379Sdas 266176379Sdas /* 267176379Sdas * XXX missing: 268176379Sdas * - tests for ld128 269176379Sdas * - tests for other rounding modes (probably won't pass for now) 270176379Sdas * - tests for large numbers that get reduced to hi+lo with lo!=0 271176379Sdas */ 272176379Sdas} 273176379Sdas 274176379Sdasint 275176379Sdasmain(int argc, char *argv[]) 276176379Sdas{ 277176379Sdas 278176379Sdas printf("1..3\n"); 279176379Sdas 280176379Sdas run_special_tests(); 281176379Sdas printf("ok 1 - trig\n"); 282176379Sdas 283176379Sdas#ifndef __i386__ 284176379Sdas run_reduction_tests(); 285176379Sdas#endif 286176379Sdas printf("ok 2 - trig\n"); 287176379Sdas 288176379Sdas#ifndef __i386__ 289176379Sdas run_accuracy_tests(); 290176379Sdas#endif 291176379Sdas printf("ok 3 - trig\n"); 292176379Sdas 293176379Sdas return (0); 294176379Sdas} 295