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