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