tools/regression/lib/msun/test-trig.c

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