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