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 <assert.h> 37 #include <complex.h> 38 #include <float.h> 39 #include <math.h> 40 #include <stdio.h> 41 42 #include "test-utils.h" 43 44 /* 45 * This is a test hook that can point to csqrtl(), _csqrt(), or to _csqrtf(). 46 * The latter two convert to float or double, respectively, and test csqrtf() 47 * and csqrt() with the same arguments. 48 */ 49 static long double complex (*t_csqrt)(long double complex); 50 51 static long double complex 52 _csqrtf(long double complex d) 53 { 54 55 return (csqrtf((float complex)d)); 56 } 57 58 static long double complex 59 _csqrt(long double complex d) 60 { 61 62 return (csqrt((double complex)d)); 63 } 64 65 #pragma STDC CX_LIMITED_RANGE OFF 66 67 /* 68 * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0. 69 * Fail an assertion if they differ. 70 */ 71 static void 72 assert_equal(long double complex d1, long double complex d2) 73 { 74 75 assert(cfpequal(d1, d2)); 76 } 77 78 /* 79 * Test csqrt for some finite arguments where the answer is exact. 80 * (We do not test if it produces correctly rounded answers when the 81 * result is inexact, nor do we check whether it throws spurious 82 * exceptions.) 83 */ 84 static void 85 test_finite(void) 86 { 87 static const double tests[] = { 88 /* csqrt(a + bI) = x + yI */ 89 /* a b x y */ 90 0, 8, 2, 2, 91 0, -8, 2, -2, 92 4, 0, 2, 0, 93 -4, 0, 0, 2, 94 3, 4, 2, 1, 95 3, -4, 2, -1, 96 -3, 4, 1, 2, 97 -3, -4, 1, -2, 98 5, 12, 3, 2, 99 7, 24, 4, 3, 100 9, 40, 5, 4, 101 11, 60, 6, 5, 102 13, 84, 7, 6, 103 33, 56, 7, 4, 104 39, 80, 8, 5, 105 65, 72, 9, 4, 106 987, 9916, 74, 67, 107 5289, 6640, 83, 40, 108 460766389075.0, 16762287900.0, 678910, 12345 109 }; 110 /* 111 * We also test some multiples of the above arguments. This 112 * array defines which multiples we use. Note that these have 113 * to be small enough to not cause overflow for float precision 114 * with all of the constants in the above table. 115 */ 116 static const double mults[] = { 117 1, 118 2, 119 3, 120 13, 121 16, 122 0x1.p30, 123 0x1.p-30, 124 }; 125 126 double a, b; 127 double x, y; 128 unsigned i, j; 129 130 for (i = 0; i < nitems(tests); i += 4) { 131 for (j = 0; j < nitems(mults); j++) { 132 a = tests[i] * mults[j] * mults[j]; 133 b = tests[i + 1] * mults[j] * mults[j]; 134 x = tests[i + 2] * mults[j]; 135 y = tests[i + 3] * mults[j]; 136 assert(t_csqrt(CMPLXL(a, b)) == CMPLXL(x, y)); 137 } 138 } 139 140 } 141 142 /* 143 * Test the handling of +/- 0. 144 */ 145 static void 146 test_zeros(void) 147 { 148 149 assert_equal(t_csqrt(CMPLXL(0.0, 0.0)), CMPLXL(0.0, 0.0)); 150 assert_equal(t_csqrt(CMPLXL(-0.0, 0.0)), CMPLXL(0.0, 0.0)); 151 assert_equal(t_csqrt(CMPLXL(0.0, -0.0)), CMPLXL(0.0, -0.0)); 152 assert_equal(t_csqrt(CMPLXL(-0.0, -0.0)), CMPLXL(0.0, -0.0)); 153 } 154 155 /* 156 * Test the handling of infinities when the other argument is not NaN. 157 */ 158 static void 159 test_infinities(void) 160 { 161 static const double vals[] = { 162 0.0, 163 -0.0, 164 42.0, 165 -42.0, 166 INFINITY, 167 -INFINITY, 168 }; 169 170 unsigned i; 171 172 for (i = 0; i < nitems(vals); i++) { 173 if (isfinite(vals[i])) { 174 assert_equal(t_csqrt(CMPLXL(-INFINITY, vals[i])), 175 CMPLXL(0.0, copysignl(INFINITY, vals[i]))); 176 assert_equal(t_csqrt(CMPLXL(INFINITY, vals[i])), 177 CMPLXL(INFINITY, copysignl(0.0, vals[i]))); 178 } 179 assert_equal(t_csqrt(CMPLXL(vals[i], INFINITY)), 180 CMPLXL(INFINITY, INFINITY)); 181 assert_equal(t_csqrt(CMPLXL(vals[i], -INFINITY)), 182 CMPLXL(INFINITY, -INFINITY)); 183 } 184 } 185 186 /* 187 * Test the handling of NaNs. 188 */ 189 static void 190 test_nans(void) 191 { 192 193 assert(creall(t_csqrt(CMPLXL(INFINITY, NAN))) == INFINITY); 194 assert(isnan(cimagl(t_csqrt(CMPLXL(INFINITY, NAN))))); 195 196 assert(isnan(creall(t_csqrt(CMPLXL(-INFINITY, NAN))))); 197 assert(isinf(cimagl(t_csqrt(CMPLXL(-INFINITY, NAN))))); 198 199 assert_equal(t_csqrt(CMPLXL(NAN, INFINITY)), 200 CMPLXL(INFINITY, INFINITY)); 201 assert_equal(t_csqrt(CMPLXL(NAN, -INFINITY)), 202 CMPLXL(INFINITY, -INFINITY)); 203 204 assert_equal(t_csqrt(CMPLXL(0.0, NAN)), CMPLXL(NAN, NAN)); 205 assert_equal(t_csqrt(CMPLXL(-0.0, NAN)), CMPLXL(NAN, NAN)); 206 assert_equal(t_csqrt(CMPLXL(42.0, NAN)), CMPLXL(NAN, NAN)); 207 assert_equal(t_csqrt(CMPLXL(-42.0, NAN)), CMPLXL(NAN, NAN)); 208 assert_equal(t_csqrt(CMPLXL(NAN, 0.0)), CMPLXL(NAN, NAN)); 209 assert_equal(t_csqrt(CMPLXL(NAN, -0.0)), CMPLXL(NAN, NAN)); 210 assert_equal(t_csqrt(CMPLXL(NAN, 42.0)), CMPLXL(NAN, NAN)); 211 assert_equal(t_csqrt(CMPLXL(NAN, -42.0)), CMPLXL(NAN, NAN)); 212 assert_equal(t_csqrt(CMPLXL(NAN, NAN)), CMPLXL(NAN, NAN)); 213 } 214 215 /* 216 * Test whether csqrt(a + bi) works for inputs that are large enough to 217 * cause overflow in hypot(a, b) + a. In this case we are using 218 * csqrt(115 + 252*I) == 14 + 9*I 219 * scaled up to near MAX_EXP. 220 */ 221 static void 222 test_overflow(int maxexp) 223 { 224 long double a, b; 225 long double complex result; 226 227 a = ldexpl(115 * 0x1p-8, maxexp); 228 b = ldexpl(252 * 0x1p-8, maxexp); 229 result = t_csqrt(CMPLXL(a, b)); 230 assert(creall(result) == ldexpl(14 * 0x1p-4, maxexp / 2)); 231 assert(cimagl(result) == ldexpl(9 * 0x1p-4, maxexp / 2)); 232 } 233 234 int 235 main(void) 236 { 237 238 printf("1..15\n"); 239 240 /* Test csqrt() */ 241 t_csqrt = _csqrt; 242 243 test_finite(); 244 printf("ok 1 - csqrt\n"); 245 246 test_zeros(); 247 printf("ok 2 - csqrt\n"); 248 249 test_infinities(); 250 printf("ok 3 - csqrt\n"); 251 252 test_nans(); 253 printf("ok 4 - csqrt\n"); 254 255 test_overflow(DBL_MAX_EXP); 256 printf("ok 5 - csqrt\n"); 257 258 /* Now test csqrtf() */ 259 t_csqrt = _csqrtf; 260 261 test_finite(); 262 printf("ok 6 - csqrt\n"); 263 264 test_zeros(); 265 printf("ok 7 - csqrt\n"); 266 267 test_infinities(); 268 printf("ok 8 - csqrt\n"); 269 270 test_nans(); 271 printf("ok 9 - csqrt\n"); 272 273 test_overflow(FLT_MAX_EXP); 274 printf("ok 10 - csqrt\n"); 275 276 /* Now test csqrtl() */ 277 t_csqrt = csqrtl; 278 279 test_finite(); 280 printf("ok 11 - csqrt\n"); 281 282 test_zeros(); 283 printf("ok 12 - csqrt\n"); 284 285 test_infinities(); 286 printf("ok 13 - csqrt\n"); 287 288 test_nans(); 289 printf("ok 14 - csqrt\n"); 290 291 test_overflow(LDBL_MAX_EXP); 292 printf("ok 15 - csqrt\n"); 293 294 return (0); 295 } 296