tools/regression/lib/msun/test-csqrt.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 <assert.h>
  35 #include <complex.h>
  36 #include <float.h>
  37 #include <math.h>
  38 #include <stdio.h>
  39
  40 #define N(i)    (sizeof(i) / sizeof((i)[0]))
  41
  42 /*
  43  * This is a test hook that can point to csqrtl(), _csqrt(), or to _csqrtf().
  44  * The latter two convert to float or double, respectively, and test csqrtf()
  45  * and csqrt() with the same arguments.
  46  */
  47 long double complex (*t_csqrt)(long double complex);
  48
  49 static long double complex
  50 _csqrtf(long double complex d)
  51 {
  52
  53         return (csqrtf((float complex)d));
  54 }
  55
  56 static long double complex
  57 _csqrt(long double complex d)
  58 {
  59
  60         return (csqrt((double complex)d));
  61 }
  62
  63 #pragma STDC CX_LIMITED_RANGE   off
  64
  65 /*
  66  * XXX gcc implements complex multiplication incorrectly. In
  67  * particular, it implements it as if the CX_LIMITED_RANGE pragma
  68  * were ON. Consequently, we need this function to form numbers
  69  * such as x + INFINITY * I, since gcc evalutes INFINITY * I as
  70  * NaN + INFINITY * I.
  71  */
  72 static inline long double complex
  73 cpackl(long double x, long double y)
  74 {
  75         long double complex z;
  76
  77         __real__ z = x;
  78         __imag__ z = y;
  79         return (z);
  80 }
  81
  82 /*
  83  * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0.
  84  * Fail an assertion if they differ.
  85  */
  86 static void
  87 assert_equal(long double complex d1, long double complex d2)
  88 {
  89
  90         if (isnan(creall(d1))) {
  91                 assert(isnan(creall(d2)));
  92         } else {
  93                 assert(creall(d1) == creall(d2));
  94                 assert(copysignl(1.0, creall(d1)) ==
  95                        copysignl(1.0, creall(d2)));
  96         }
  97         if (isnan(cimagl(d1))) {
  98                 assert(isnan(cimagl(d2)));
  99         } else {
 100                 assert(cimagl(d1) == cimagl(d2));
 101                 assert(copysignl(1.0, cimagl(d1)) ==
 102                        copysignl(1.0, cimagl(d2)));
 103         }
 104 }
 105
 106 /*
 107  * Test csqrt for some finite arguments where the answer is exact.
 108  * (We do not test if it produces correctly rounded answers when the
 109  * result is inexact, nor do we check whether it throws spurious
 110  * exceptions.)
 111  */
 112 static void
 113 test_finite()
 114 {
 115         static const double tests[] = {
 116              /* csqrt(a + bI) = x + yI */
 117              /* a       b       x       y */
 118                 0,      8,      2,      2,
 119                 0,      -8,     2,      -2,
 120                 4,      0,      2,      0,
 121                 -4,     0,      0,      2,
 122                 3,      4,      2,      1,
 123                 3,      -4,     2,      -1,
 124                 -3,     4,      1,      2,
 125                 -3,     -4,     1,      -2,
 126                 5,      12,     3,      2,
 127                 7,      24,     4,      3,
 128                 9,      40,     5,      4,
 129                 11,     60,     6,      5,
 130                 13,     84,     7,      6,
 131                 33,     56,     7,      4,
 132                 39,     80,     8,      5,
 133                 65,     72,     9,      4,
 134                 987,    9916,   74,     67,
 135                 5289,   6640,   83,     40,
 136                 460766389075.0, 16762287900.0, 678910, 12345
 137         };
 138         /*
 139          * We also test some multiples of the above arguments. This
 140          * array defines which multiples we use. Note that these have
 141          * to be small enough to not cause overflow for float precision
 142          * with all of the constants in the above table.
 143          */
 144         static const double mults[] = {
 145                 1,
 146                 2,
 147                 3,
 148                 13,
 149                 16,
 150                 0x1.p30,
 151                 0x1.p-30,
 152         };
 153
 154         double a, b;
 155         double x, y;
 156         int i, j;
 157
 158         for (i = 0; i < N(tests); i += 4) {
 159                 for (j = 0; j < N(mults); j++) {
 160                         a = tests[i] * mults[j] * mults[j];
 161                         b = tests[i + 1] * mults[j] * mults[j];
 162                         x = tests[i + 2] * mults[j];
 163                         y = tests[i + 3] * mults[j];
 164                         assert(t_csqrt(cpackl(a, b)) == cpackl(x, y));
 165                 }
 166         }
 167
 168 }
 169
 170 /*
 171  * Test the handling of +/- 0.
 172  */
 173 static void
 174 test_zeros()
 175 {
 176
 177         assert_equal(t_csqrt(cpackl(0.0, 0.0)), cpackl(0.0, 0.0));
 178         assert_equal(t_csqrt(cpackl(-0.0, 0.0)), cpackl(0.0, 0.0));
 179         assert_equal(t_csqrt(cpackl(0.0, -0.0)), cpackl(0.0, -0.0));
 180         assert_equal(t_csqrt(cpackl(-0.0, -0.0)), cpackl(0.0, -0.0));
 181 }
 182
 183 /*
 184  * Test the handling of infinities when the other argument is not NaN.
 185  */
 186 static void
 187 test_infinities()
 188 {
 189         static const double vals[] = {
 190                 0.0,
 191                 -0.0,
 192                 42.0,
 193                 -42.0,
 194                 INFINITY,
 195                 -INFINITY,
 196         };
 197
 198         int i;
 199
 200         for (i = 0; i < N(vals); i++) {
 201                 if (isfinite(vals[i])) {
 202                         assert_equal(t_csqrt(cpackl(-INFINITY, vals[i])),
 203                             cpackl(0.0, copysignl(INFINITY, vals[i])));
 204                         assert_equal(t_csqrt(cpackl(INFINITY, vals[i])),
 205                             cpackl(INFINITY, copysignl(0.0, vals[i])));
 206                 }
 207                 assert_equal(t_csqrt(cpackl(vals[i], INFINITY)),
 208                     cpackl(INFINITY, INFINITY));
 209                 assert_equal(t_csqrt(cpackl(vals[i], -INFINITY)),
 210                     cpackl(INFINITY, -INFINITY));
 211         }
 212 }
 213
 214 /*
 215  * Test the handling of NaNs.
 216  */
 217 static void
 218 test_nans()
 219 {
 220
 221         assert(creall(t_csqrt(cpackl(INFINITY, NAN))) == INFINITY);
 222         assert(isnan(cimagl(t_csqrt(cpackl(INFINITY, NAN)))));
 223
 224         assert(isnan(creall(t_csqrt(cpackl(-INFINITY, NAN)))));
 225         assert(isinf(cimagl(t_csqrt(cpackl(-INFINITY, NAN)))));
 226
 227         assert_equal(t_csqrt(cpackl(NAN, INFINITY)),
 228                      cpackl(INFINITY, INFINITY));
 229         assert_equal(t_csqrt(cpackl(NAN, -INFINITY)),
 230                      cpackl(INFINITY, -INFINITY));
 231
 232         assert_equal(t_csqrt(cpackl(0.0, NAN)), cpackl(NAN, NAN));
 233         assert_equal(t_csqrt(cpackl(-0.0, NAN)), cpackl(NAN, NAN));
 234         assert_equal(t_csqrt(cpackl(42.0, NAN)), cpackl(NAN, NAN));
 235         assert_equal(t_csqrt(cpackl(-42.0, NAN)), cpackl(NAN, NAN));
 236         assert_equal(t_csqrt(cpackl(NAN, 0.0)), cpackl(NAN, NAN));
 237         assert_equal(t_csqrt(cpackl(NAN, -0.0)), cpackl(NAN, NAN));
 238         assert_equal(t_csqrt(cpackl(NAN, 42.0)), cpackl(NAN, NAN));
 239         assert_equal(t_csqrt(cpackl(NAN, -42.0)), cpackl(NAN, NAN));
 240         assert_equal(t_csqrt(cpackl(NAN, NAN)), cpackl(NAN, NAN));
 241 }
 242
 243 /*
 244  * Test whether csqrt(a + bi) works for inputs that are large enough to
 245  * cause overflow in hypot(a, b) + a. In this case we are using
 246  *      csqrt(115 + 252*I) == 14 + 9*I
 247  * scaled up to near MAX_EXP.
 248  */
 249 static void
 250 test_overflow(int maxexp)
 251 {
 252         long double a, b;
 253         long double complex result;
 254
 255         a = ldexpl(115 * 0x1p-8, maxexp);
 256         b = ldexpl(252 * 0x1p-8, maxexp);
 257         result = t_csqrt(cpackl(a, b));
 258         assert(creall(result) == ldexpl(14 * 0x1p-4, maxexp / 2));
 259         assert(cimagl(result) == ldexpl(9 * 0x1p-4, maxexp / 2));
 260 }
 261
 262 int
 263 main(int argc, char *argv[])
 264 {
 265
 266         printf("1..15\n");
 267
 268         /* Test csqrt() */
 269         t_csqrt = _csqrt;
 270
 271         test_finite();
 272         printf("ok 1 - csqrt\n");
 273
 274         test_zeros();
 275         printf("ok 2 - csqrt\n");
 276
 277         test_infinities();
 278         printf("ok 3 - csqrt\n");
 279
 280         test_nans();
 281         printf("ok 4 - csqrt\n");
 282
 283         test_overflow(DBL_MAX_EXP);
 284         printf("ok 5 - csqrt\n");
 285
 286         /* Now test csqrtf() */
 287         t_csqrt = _csqrtf;
 288
 289         test_finite();
 290         printf("ok 6 - csqrt\n");
 291
 292         test_zeros();
 293         printf("ok 7 - csqrt\n");
 294
 295         test_infinities();
 296         printf("ok 8 - csqrt\n");
 297
 298         test_nans();
 299         printf("ok 9 - csqrt\n");
 300
 301         test_overflow(FLT_MAX_EXP);
 302         printf("ok 10 - csqrt\n");
 303
 304         /* Now test csqrtl() */
 305         t_csqrt = csqrtl;
 306
 307         test_finite();
 308         printf("ok 11 - csqrt\n");
 309
 310         test_zeros();
 311         printf("ok 12 - csqrt\n");
 312
 313         test_infinities();
 314         printf("ok 13 - csqrt\n");
 315
 316         test_nans();
 317         printf("ok 14 - csqrt\n");
 318
 319         test_overflow(LDBL_MAX_EXP);
 320         printf("ok 15 - csqrt\n");
 321
 322         return (0);
 323 }