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