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