1 /* $OpenBSD: csqrt_test.c,v 1.1 2021/10/22 18:00:22 mbuhl Exp $ */ 2 /*- 3 * Copyright (c) 2007 David Schultz <das@FreeBSD.org> 4 * All rights reserved. 5 * 6 * Redistribution and use in source and binary forms, with or without 7 * modification, are permitted provided that the following conditions 8 * are met: 9 * 1. Redistributions of source code must retain the above copyright 10 * notice, this list of conditions and the following disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 18 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 19 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 20 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 21 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 25 * SUCH DAMAGE. 26 */ 27 28 #include "macros.h" 29 30 #include <sys/cdefs.h> 31 __FBSDID("$FreeBSD$"); 32 33 #include <sys/param.h> 34 35 #include <complex.h> 36 #include <float.h> 37 #include <math.h> 38 #include <stdio.h> 39 40 #include "test-utils.h" 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 static 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 * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0. 67 * Fail an assertion if they differ. 68 */ 69 #define assert_equal(d1, d2) CHECK_CFPEQUAL_CS(d1, d2, CS_BOTH) 70 71 /* 72 * Test csqrt for some finite arguments where the answer is exact. 73 * (We do not test if it produces correctly rounded answers when the 74 * result is inexact, nor do we check whether it throws spurious 75 * exceptions.) 76 */ 77 static void 78 test_finite(void) 79 { 80 static const double tests[] = { 81 /* csqrt(a + bI) = x + yI */ 82 /* a b x y */ 83 0, 8, 2, 2, 84 0, -8, 2, -2, 85 4, 0, 2, 0, 86 -4, 0, 0, 2, 87 3, 4, 2, 1, 88 3, -4, 2, -1, 89 -3, 4, 1, 2, 90 -3, -4, 1, -2, 91 5, 12, 3, 2, 92 7, 24, 4, 3, 93 9, 40, 5, 4, 94 11, 60, 6, 5, 95 13, 84, 7, 6, 96 33, 56, 7, 4, 97 39, 80, 8, 5, 98 65, 72, 9, 4, 99 987, 9916, 74, 67, 100 5289, 6640, 83, 40, 101 460766389075.0, 16762287900.0, 678910, 12345 102 }; 103 /* 104 * We also test some multiples of the above arguments. This 105 * array defines which multiples we use. Note that these have 106 * to be small enough to not cause overflow for float precision 107 * with all of the constants in the above table. 108 */ 109 static const double mults[] = { 110 1, 111 2, 112 3, 113 13, 114 16, 115 0x1.p30, 116 0x1.p-30, 117 }; 118 119 double a, b; 120 double x, y; 121 unsigned i, j; 122 123 for (i = 0; i < nitems(tests); i += 4) { 124 for (j = 0; j < nitems(mults); j++) { 125 a = tests[i] * mults[j] * mults[j]; 126 b = tests[i + 1] * mults[j] * mults[j]; 127 x = tests[i + 2] * mults[j]; 128 y = tests[i + 3] * mults[j]; 129 ATF_CHECK(t_csqrt(CMPLXL(a, b)) == CMPLXL(x, y)); 130 } 131 } 132 133 } 134 135 /* 136 * Test the handling of +/- 0. 137 */ 138 static void 139 test_zeros(void) 140 { 141 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 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 } 147 148 /* 149 * Test the handling of infinities when the other argument is not NaN. 150 */ 151 static void 152 test_infinities(void) 153 { 154 static const double vals[] = { 155 0.0, 156 -0.0, 157 42.0, 158 -42.0, 159 INFINITY, 160 -INFINITY, 161 }; 162 163 unsigned i; 164 165 for (i = 0; i < nitems(vals); i++) { 166 if (isfinite(vals[i])) { 167 assert_equal(t_csqrt(CMPLXL(-INFINITY, vals[i])), 168 CMPLXL(0.0, copysignl(INFINITY, vals[i]))); 169 assert_equal(t_csqrt(CMPLXL(INFINITY, vals[i])), 170 CMPLXL(INFINITY, copysignl(0.0, vals[i]))); 171 } 172 assert_equal(t_csqrt(CMPLXL(vals[i], INFINITY)), 173 CMPLXL(INFINITY, INFINITY)); 174 assert_equal(t_csqrt(CMPLXL(vals[i], -INFINITY)), 175 CMPLXL(INFINITY, -INFINITY)); 176 } 177 } 178 179 /* 180 * Test the handling of NaNs. 181 */ 182 static void 183 test_nans(void) 184 { 185 186 ATF_CHECK(creall(t_csqrt(CMPLXL(INFINITY, NAN))) == INFINITY); 187 ATF_CHECK(isnan(cimagl(t_csqrt(CMPLXL(INFINITY, NAN))))); 188 189 ATF_CHECK(isnan(creall(t_csqrt(CMPLXL(-INFINITY, NAN))))); 190 ATF_CHECK(isinf(cimagl(t_csqrt(CMPLXL(-INFINITY, NAN))))); 191 192 assert_equal(t_csqrt(CMPLXL(NAN, INFINITY)), 193 CMPLXL(INFINITY, INFINITY)); 194 assert_equal(t_csqrt(CMPLXL(NAN, -INFINITY)), 195 CMPLXL(INFINITY, -INFINITY)); 196 197 assert_equal(t_csqrt(CMPLXL(0.0, NAN)), CMPLXL(NAN, NAN)); 198 assert_equal(t_csqrt(CMPLXL(-0.0, NAN)), CMPLXL(NAN, NAN)); 199 assert_equal(t_csqrt(CMPLXL(42.0, NAN)), CMPLXL(NAN, NAN)); 200 assert_equal(t_csqrt(CMPLXL(-42.0, NAN)), CMPLXL(NAN, NAN)); 201 assert_equal(t_csqrt(CMPLXL(NAN, 0.0)), CMPLXL(NAN, NAN)); 202 assert_equal(t_csqrt(CMPLXL(NAN, -0.0)), CMPLXL(NAN, NAN)); 203 assert_equal(t_csqrt(CMPLXL(NAN, 42.0)), CMPLXL(NAN, NAN)); 204 assert_equal(t_csqrt(CMPLXL(NAN, -42.0)), CMPLXL(NAN, NAN)); 205 assert_equal(t_csqrt(CMPLXL(NAN, NAN)), CMPLXL(NAN, NAN)); 206 } 207 208 /* 209 * Test whether csqrt(a + bi) works for inputs that are large enough to 210 * cause overflow in hypot(a, b) + a. Each of the tests is scaled up to 211 * near MAX_EXP. 212 */ 213 static void 214 test_overflow(int maxexp) 215 { 216 long double a, b; 217 long double complex result; 218 int exp, i; 219 220 ATF_CHECK(maxexp > 0 && maxexp % 2 == 0); 221 222 for (i = 0; i < 4; i++) { 223 exp = maxexp - 2 * i; 224 225 /* csqrt(115 + 252*I) == 14 + 9*I */ 226 a = ldexpl(115 * 0x1p-8, exp); 227 b = ldexpl(252 * 0x1p-8, exp); 228 result = t_csqrt(CMPLXL(a, b)); 229 ATF_CHECK_EQ(creall(result), ldexpl(14 * 0x1p-4, exp / 2)); 230 ATF_CHECK_EQ(cimagl(result), ldexpl(9 * 0x1p-4, exp / 2)); 231 232 /* csqrt(-11 + 60*I) = 5 + 6*I */ 233 a = ldexpl(-11 * 0x1p-6, exp); 234 b = ldexpl(60 * 0x1p-6, exp); 235 result = t_csqrt(CMPLXL(a, b)); 236 ATF_CHECK_EQ(creall(result), ldexpl(5 * 0x1p-3, exp / 2)); 237 ATF_CHECK_EQ(cimagl(result), ldexpl(6 * 0x1p-3, exp / 2)); 238 239 /* csqrt(225 + 0*I) == 15 + 0*I */ 240 a = ldexpl(225 * 0x1p-8, exp); 241 b = 0; 242 result = t_csqrt(CMPLXL(a, b)); 243 ATF_CHECK_EQ(creall(result), ldexpl(15 * 0x1p-4, exp / 2)); 244 ATF_CHECK_EQ(cimagl(result), 0); 245 } 246 } 247 248 /* 249 * Test that precision is maintained for some large squares. Set all or 250 * some bits in the lower mantdig/2 bits, square the number, and try to 251 * recover the sqrt. Note: 252 * (x + xI)**2 = 2xxI 253 */ 254 static void 255 test_precision(int maxexp, int mantdig) 256 { 257 long double b, x; 258 long double complex result; 259 #if LDBL_MANT_DIG <= 64 260 typedef uint64_t ldbl_mant_type; 261 #elif LDBL_MANT_DIG <= 128 262 typedef __uint128_t ldbl_mant_type; 263 #else 264 #error "Unsupported long double format" 265 #endif 266 ldbl_mant_type mantbits, sq_mantbits; 267 int exp, i; 268 269 ATF_REQUIRE(maxexp > 0 && maxexp % 2 == 0); 270 ATF_REQUIRE(mantdig <= LDBL_MANT_DIG); 271 mantdig = rounddown(mantdig, 2); 272 273 for (exp = 0; exp <= maxexp; exp += 2) { 274 mantbits = ((ldbl_mant_type)1 << (mantdig / 2)) - 1; 275 for (i = 0; i < 100 && 276 mantbits > ((ldbl_mant_type)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 CHECK_FPEQUAL(b, x * x * 2); 293 result = t_csqrt(CMPLXL(0, b)); 294 CHECK_FPEQUAL(x, creall(result)); 295 CHECK_FPEQUAL(x, cimagl(result)); 296 } 297 } 298 } 299 300 ATF_TC_WITHOUT_HEAD(csqrt); 301 ATF_TC_BODY(csqrt, tc) 302 { 303 /* Test csqrt() */ 304 t_csqrt = _csqrt; 305 306 test_finite(); 307 308 test_zeros(); 309 310 test_infinities(); 311 312 test_nans(); 313 314 test_overflow(DBL_MAX_EXP); 315 316 test_precision(DBL_MAX_EXP, DBL_MANT_DIG); 317 } 318 319 ATF_TC_WITHOUT_HEAD(csqrtf); 320 ATF_TC_BODY(csqrtf, tc) 321 { 322 /* Now test csqrtf() */ 323 t_csqrt = _csqrtf; 324 325 test_finite(); 326 327 test_zeros(); 328 329 test_infinities(); 330 331 test_nans(); 332 333 test_overflow(FLT_MAX_EXP); 334 335 test_precision(FLT_MAX_EXP, FLT_MANT_DIG); 336 } 337 338 ATF_TC_WITHOUT_HEAD(csqrtl); 339 ATF_TC_BODY(csqrtl, tc) 340 { 341 /* Now test csqrtl() */ 342 t_csqrt = csqrtl; 343 344 test_finite(); 345 346 test_zeros(); 347 348 test_infinities(); 349 350 test_nans(); 351 352 test_overflow(LDBL_MAX_EXP); 353 354 /* i386 is configured to use 53-bit rounding precision for long double. */ 355 #ifndef __i386__ 356 test_precision(LDBL_MAX_EXP, LDBL_MANT_DIG); 357 #else 358 test_precision(LDBL_MAX_EXP, DBL_MANT_DIG); 359 #endif 360 } 361 362 ATF_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