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