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 } 373