lib/msun/tests/csqrt_test.c

   1 /*-
   2  * Copyright (c) 2007 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 csqrt{,f}()
  29  */
  30
  31 #include <sys/cdefs.h>
  32 __FBSDID("$FreeBSD$");
  33
  34 #include <sys/param.h>
  35
  36 #include <complex.h>
  37 #include <float.h>
  38 #include <math.h>
  39 #include <stdio.h>
  40
  41 #include "test-utils.h"
  42
  43 /*
  44  * This is a test hook that can point to csqrtl(), _csqrt(), or to _csqrtf().
  45  * The latter two convert to float or double, respectively, and test csqrtf()
  46  * and csqrt() with the same arguments.
  47  */
  48 static long double complex (*t_csqrt)(long double complex);
  49
  50 static long double complex
  51 _csqrtf(long double complex d)
  52 {
  53
  54         return (csqrtf((float complex)d));
  55 }
  56
  57 static long double complex
  58 _csqrt(long double complex d)
  59 {
  60
  61         return (csqrt((double complex)d));
  62 }
  63
  64 #pragma STDC CX_LIMITED_RANGE   OFF
  65
  66 /*
  67  * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0.
  68  * Fail an assertion if they differ.
  69  */
  70 #define assert_equal(d1, d2) CHECK_CFPEQUAL_CS(d1, d2, CS_BOTH)
  71
  72 /*
  73  * Test csqrt for some finite arguments where the answer is exact.
  74  * (We do not test if it produces correctly rounded answers when the
  75  * result is inexact, nor do we check whether it throws spurious
  76  * exceptions.)
  77  */
  78 static void
  79 test_finite(void)
  80 {
  81         static const double tests[] = {
  82              /* csqrt(a + bI) = x + yI */
  83              /* a       b       x       y */
  84                 0,      8,      2,      2,
  85                 0,      -8,     2,      -2,
  86                 4,      0,      2,      0,
  87                 -4,     0,      0,      2,
  88                 3,      4,      2,      1,
  89                 3,      -4,     2,      -1,
  90                 -3,     4,      1,      2,
  91                 -3,     -4,     1,      -2,
  92                 5,      12,     3,      2,
  93                 7,      24,     4,      3,
  94                 9,      40,     5,      4,
  95                 11,     60,     6,      5,
  96                 13,     84,     7,      6,
  97                 33,     56,     7,      4,
  98                 39,     80,     8,      5,
  99                 65,     72,     9,      4,
 100                 987,    9916,   74,     67,
 101                 5289,   6640,   83,     40,
 102                 460766389075.0, 16762287900.0, 678910, 12345
 103         };
 104         /*
 105          * We also test some multiples of the above arguments. This
 106          * array defines which multiples we use. Note that these have
 107          * to be small enough to not cause overflow for float precision
 108          * with all of the constants in the above table.
 109          */
 110         static const double mults[] = {
 111                 1,
 112                 2,
 113                 3,
 114                 13,
 115                 16,
 116                 0x1.p30,
 117                 0x1.p-30,
 118         };
 119
 120         double a, b;
 121         double x, y;
 122         unsigned i, j;
 123
 124         for (i = 0; i < nitems(tests); i += 4) {
 125                 for (j = 0; j < nitems(mults); j++) {
 126                         a = tests[i] * mults[j] * mults[j];
 127                         b = tests[i + 1] * mults[j] * mults[j];
 128                         x = tests[i + 2] * mults[j];
 129                         y = tests[i + 3] * mults[j];
 130                         ATF_CHECK(t_csqrt(CMPLXL(a, b)) == CMPLXL(x, y));
 131                 }
 132         }
 133
 134 }
 135
 136 /*
 137  * Test the handling of +/- 0.
 138  */
 139 static void
 140 test_zeros(void)
 141 {
 142
 143         assert_equal(t_csqrt(CMPLXL(0.0, 0.0)), CMPLXL(0.0, 0.0));
 144         assert_equal(t_csqrt(CMPLXL(-0.0, 0.0)), CMPLXL(0.0, 0.0));
 145         assert_equal(t_csqrt(CMPLXL(0.0, -0.0)), CMPLXL(0.0, -0.0));
 146         assert_equal(t_csqrt(CMPLXL(-0.0, -0.0)), CMPLXL(0.0, -0.0));
 147 }
 148
 149 /*
 150  * Test the handling of infinities when the other argument is not NaN.
 151  */
 152 static void
 153 test_infinities(void)
 154 {
 155         static const double vals[] = {
 156                 0.0,
 157                 -0.0,
 158                 42.0,
 159                 -42.0,
 160                 INFINITY,
 161                 -INFINITY,
 162         };
 163
 164         unsigned i;
 165
 166         for (i = 0; i < nitems(vals); i++) {
 167                 if (isfinite(vals[i])) {
 168                         assert_equal(t_csqrt(CMPLXL(-INFINITY, vals[i])),
 169                             CMPLXL(0.0, copysignl(INFINITY, vals[i])));
 170                         assert_equal(t_csqrt(CMPLXL(INFINITY, vals[i])),
 171                             CMPLXL(INFINITY, copysignl(0.0, vals[i])));
 172                 }
 173                 assert_equal(t_csqrt(CMPLXL(vals[i], INFINITY)),
 174                     CMPLXL(INFINITY, INFINITY));
 175                 assert_equal(t_csqrt(CMPLXL(vals[i], -INFINITY)),
 176                     CMPLXL(INFINITY, -INFINITY));
 177         }
 178 }
 179
 180 /*
 181  * Test the handling of NaNs.
 182  */
 183 static void
 184 test_nans(void)
 185 {
 186
 187         ATF_CHECK(creall(t_csqrt(CMPLXL(INFINITY, NAN))) == INFINITY);
 188         ATF_CHECK(isnan(cimagl(t_csqrt(CMPLXL(INFINITY, NAN)))));
 189
 190         ATF_CHECK(isnan(creall(t_csqrt(CMPLXL(-INFINITY, NAN)))));
 191         ATF_CHECK(isinf(cimagl(t_csqrt(CMPLXL(-INFINITY, NAN)))));
 192
 193         assert_equal(t_csqrt(CMPLXL(NAN, INFINITY)),
 194                      CMPLXL(INFINITY, INFINITY));
 195         assert_equal(t_csqrt(CMPLXL(NAN, -INFINITY)),
 196                      CMPLXL(INFINITY, -INFINITY));
 197
 198         assert_equal(t_csqrt(CMPLXL(0.0, NAN)), CMPLXL(NAN, NAN));
 199         assert_equal(t_csqrt(CMPLXL(-0.0, NAN)), CMPLXL(NAN, NAN));
 200         assert_equal(t_csqrt(CMPLXL(42.0, NAN)), CMPLXL(NAN, NAN));
 201         assert_equal(t_csqrt(CMPLXL(-42.0, NAN)), CMPLXL(NAN, NAN));
 202         assert_equal(t_csqrt(CMPLXL(NAN, 0.0)), CMPLXL(NAN, NAN));
 203         assert_equal(t_csqrt(CMPLXL(NAN, -0.0)), CMPLXL(NAN, NAN));
 204         assert_equal(t_csqrt(CMPLXL(NAN, 42.0)), CMPLXL(NAN, NAN));
 205         assert_equal(t_csqrt(CMPLXL(NAN, -42.0)), CMPLXL(NAN, NAN));
 206         assert_equal(t_csqrt(CMPLXL(NAN, NAN)), CMPLXL(NAN, NAN));
 207 }
 208
 209 /*
 210  * Test whether csqrt(a + bi) works for inputs that are large enough to
 211  * cause overflow in hypot(a, b) + a.  Each of the tests is scaled up to
 212  * near MAX_EXP.
 213  */
 214 static void
 215 test_overflow(int maxexp)
 216 {
 217         long double a, b;
 218         long double complex result;
 219         int exp, i;
 220
 221         ATF_CHECK(maxexp > 0 && maxexp % 2 == 0);
 222
 223         for (i = 0; i < 4; i++) {
 224                 exp = maxexp - 2 * i;
 225
 226                 /* csqrt(115 + 252*I) == 14 + 9*I */
 227                 a = ldexpl(115 * 0x1p-8, exp);
 228                 b = ldexpl(252 * 0x1p-8, exp);
 229                 result = t_csqrt(CMPLXL(a, b));
 230                 ATF_CHECK_EQ(creall(result), ldexpl(14 * 0x1p-4, exp / 2));
 231                 ATF_CHECK_EQ(cimagl(result), ldexpl(9 * 0x1p-4, exp / 2));
 232
 233                 /* csqrt(-11 + 60*I) = 5 + 6*I */
 234                 a = ldexpl(-11 * 0x1p-6, exp);
 235                 b = ldexpl(60 * 0x1p-6, exp);
 236                 result = t_csqrt(CMPLXL(a, b));
 237                 ATF_CHECK_EQ(creall(result), ldexpl(5 * 0x1p-3, exp / 2));
 238                 ATF_CHECK_EQ(cimagl(result), ldexpl(6 * 0x1p-3, exp / 2));
 239
 240                 /* csqrt(225 + 0*I) == 15 + 0*I */
 241                 a = ldexpl(225 * 0x1p-8, exp);
 242                 b = 0;
 243                 result = t_csqrt(CMPLXL(a, b));
 244                 ATF_CHECK_EQ(creall(result), ldexpl(15 * 0x1p-4, exp / 2));
 245                 ATF_CHECK_EQ(cimagl(result), 0);
 246         }
 247 }
 248
 249 /*
 250  * Test that precision is maintained for some large squares.  Set all or
 251  * some bits in the lower mantdig/2 bits, square the number, and try to
 252  * recover the sqrt.  Note:
 253  *      (x + xI)**2 = 2xxI
 254  */
 255 static void
 256 test_precision(int maxexp, int mantdig)
 257 {
 258         long double b, x;
 259         long double complex result;
 260 #if LDBL_MANT_DIG <= 64
 261         typedef uint64_t ldbl_mant_type;
 262 #elif LDBL_MANT_DIG <= 128
 263         typedef __uint128_t ldbl_mant_type;
 264 #else
 265 #error "Unsupported long double format"
 266 #endif
 267         ldbl_mant_type mantbits, sq_mantbits;
 268         int exp, i;
 269
 270         ATF_REQUIRE(maxexp > 0 && maxexp % 2 == 0);
 271         ATF_REQUIRE(mantdig <= LDBL_MANT_DIG);
 272         mantdig = rounddown(mantdig, 2);
 273
 274         for (exp = 0; exp <= maxexp; exp += 2) {
 275                 mantbits = ((ldbl_mant_type)1 << (mantdig / 2)) - 1;
 276                 for (i = 0; i < 100 &&
 277                      mantbits > ((ldbl_mant_type)1 << (mantdig / 2 - 1));
 278                      i++, mantbits--) {
 279                         sq_mantbits = mantbits * mantbits;
 280                         /*
 281                          * sq_mantibts is a mantdig-bit number.  Divide by
 282                          * 2**mantdig to normalize it to [0.5, 1), where,
 283                          * note, the binary power will be -1.  Raise it by
 284                          * 2**exp for the test.  exp is even.  Lower it by
 285                          * one to reach a final binary power which is also
 286                          * even.  The result should be exactly
 287                          * representable, given that mantdig is less than or
 288                          * equal to the available precision.
 289                          */
 290                         b = ldexpl((long double)sq_mantbits,
 291                             exp - 1 - mantdig);
 292                         x = ldexpl(mantbits, (exp - 2 - mantdig) / 2);
 293                         CHECK_FPEQUAL(b, x * x * 2);
 294                         result = t_csqrt(CMPLXL(0, b));
 295                         CHECK_FPEQUAL(x, creall(result));
 296                         CHECK_FPEQUAL(x, cimagl(result));
 297                 }
 298         }
 299 }
 300
 301 ATF_TC_WITHOUT_HEAD(csqrt);
 302 ATF_TC_BODY(csqrt, tc)
 303 {
 304         /* Test csqrt() */
 305         t_csqrt = _csqrt;
 306
 307         test_finite();
 308
 309         test_zeros();
 310
 311         test_infinities();
 312
 313         test_nans();
 314
 315         test_overflow(DBL_MAX_EXP);
 316
 317         test_precision(DBL_MAX_EXP, DBL_MANT_DIG);
 318 }
 319
 320 ATF_TC_WITHOUT_HEAD(csqrtf);
 321 ATF_TC_BODY(csqrtf, tc)
 322 {
 323         /* Now test csqrtf() */
 324         t_csqrt = _csqrtf;
 325
 326         test_finite();
 327
 328         test_zeros();
 329
 330         test_infinities();
 331
 332         test_nans();
 333
 334         test_overflow(FLT_MAX_EXP);
 335
 336         test_precision(FLT_MAX_EXP, FLT_MANT_DIG);
 337 }
 338
 339 ATF_TC_WITHOUT_HEAD(csqrtl);
 340 ATF_TC_BODY(csqrtl, tc)
 341 {
 342         /* Now test csqrtl() */
 343         t_csqrt = csqrtl;
 344
 345         test_finite();
 346
 347         test_zeros();
 348
 349         test_infinities();
 350
 351         test_nans();
 352
 353         test_overflow(LDBL_MAX_EXP);
 354
 355         /* i386 is configured to use 53-bit rounding precision for long double. */
 356         test_precision(LDBL_MAX_EXP,
 357 #ifndef __i386__
 358             LDBL_MANT_DIG
 359 #else
 360             DBL_MANT_DIG
 361 #endif
 362             );
 363 }
 364
 365 ATF_TP_ADD_TCS(tp)
 366 {
 367         ATF_TP_ADD_TC(tp, csqrt);
 368         ATF_TP_ADD_TC(tp, csqrtf);
 369         ATF_TP_ADD_TC(tp, csqrtl);
 370
 371         return (atf_no_error());
 372 }