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