csqrt_test.c revision 174619
1174619Sdas/*- 2174619Sdas * Copyright (c) 2007 David Schultz <das@FreeBSD.org> 3174619Sdas * All rights reserved. 4174619Sdas * 5174619Sdas * Redistribution and use in source and binary forms, with or without 6174619Sdas * modification, are permitted provided that the following conditions 7174619Sdas * are met: 8174619Sdas * 1. Redistributions of source code must retain the above copyright 9174619Sdas * notice, this list of conditions and the following disclaimer. 10174619Sdas * 2. Redistributions in binary form must reproduce the above copyright 11174619Sdas * notice, this list of conditions and the following disclaimer in the 12174619Sdas * documentation and/or other materials provided with the distribution. 13174619Sdas * 14174619Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15174619Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16174619Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17174619Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18174619Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19174619Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20174619Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21174619Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22174619Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23174619Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24174619Sdas * SUCH DAMAGE. 25174619Sdas */ 26174619Sdas 27174619Sdas/* 28174619Sdas * Tests for csqrt{,f}() 29174619Sdas */ 30174619Sdas 31174619Sdas#include <sys/cdefs.h> 32174619Sdas__FBSDID("$FreeBSD: head/tools/regression/lib/msun/test-csqrt.c 174619 2007-12-15 09:16:26Z das $"); 33174619Sdas 34174619Sdas#include <assert.h> 35174619Sdas#include <complex.h> 36174619Sdas#include <float.h> 37174619Sdas#include <math.h> 38174619Sdas#include <stdio.h> 39174619Sdas 40174619Sdas#define N(i) (sizeof(i) / sizeof((i)[0])) 41174619Sdas 42174619Sdas/* 43174619Sdas * This is a test hook that can point to csqrt(), or to _csqrtf(), 44174619Sdas * which converts to float and tests csqrtf() with the same arguments. 45174619Sdas */ 46174619Sdasdouble complex (*t_csqrt)(double complex); 47174619Sdas 48174619Sdasstatic double complex 49174619Sdas_csqrtf(double complex d) 50174619Sdas{ 51174619Sdas 52174619Sdas return (csqrtf((float complex)d)); 53174619Sdas} 54174619Sdas 55174619Sdas#pragma STDC CX_LIMITED_RANGE off 56174619Sdas 57174619Sdas/* 58174619Sdas * XXX gcc implements complex multiplication incorrectly. In 59174619Sdas * particular, it implements it as if the CX_LIMITED_RANGE pragma 60174619Sdas * were ON. Consequently, we need this function to form numbers 61174619Sdas * such as x + INFINITY * I, since gcc evalutes INFINITY * I as 62174619Sdas * NaN + INFINITY * I. 63174619Sdas */ 64174619Sdasstatic inline double complex 65174619Sdascpack(double x, double y) 66174619Sdas{ 67174619Sdas double complex z; 68174619Sdas 69174619Sdas __real__ z = x; 70174619Sdas __imag__ z = y; 71174619Sdas return (z); 72174619Sdas} 73174619Sdas 74174619Sdas/* 75174619Sdas * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0. 76174619Sdas * Fail an assertion if they differ. 77174619Sdas */ 78174619Sdasstatic void 79174619Sdasassert_equal(double complex d1, double complex d2) 80174619Sdas{ 81174619Sdas 82174619Sdas if (isnan(creal(d1))) { 83174619Sdas assert(isnan(creal(d2))); 84174619Sdas } else { 85174619Sdas assert(creal(d1) == creal(d2)); 86174619Sdas assert(copysign(1.0, creal(d1)) == copysign(1.0, creal(d2))); 87174619Sdas } 88174619Sdas if (isnan(cimag(d1))) { 89174619Sdas assert(isnan(cimag(d2))); 90174619Sdas } else { 91174619Sdas assert(cimag(d1) == cimag(d2)); 92174619Sdas assert(copysign(1.0, cimag(d1)) == copysign(1.0, cimag(d2))); 93174619Sdas } 94174619Sdas} 95174619Sdas 96174619Sdas/* 97174619Sdas * Test csqrt for some finite arguments where the answer is exact. 98174619Sdas * (We do not test if it produces correctly rounded answers when the 99174619Sdas * result is inexact, nor do we check whether it throws spurious 100174619Sdas * exceptions.) 101174619Sdas */ 102174619Sdasstatic void 103174619Sdastest_finite() 104174619Sdas{ 105174619Sdas static const double tests[] = { 106174619Sdas /* csqrt(a + bI) = x + yI */ 107174619Sdas /* a b x y */ 108174619Sdas 0, 8, 2, 2, 109174619Sdas 0, -8, 2, -2, 110174619Sdas 4, 0, 2, 0, 111174619Sdas -4, 0, 0, 2, 112174619Sdas 3, 4, 2, 1, 113174619Sdas 3, -4, 2, -1, 114174619Sdas -3, 4, 1, 2, 115174619Sdas -3, -4, 1, -2, 116174619Sdas 5, 12, 3, 2, 117174619Sdas 7, 24, 4, 3, 118174619Sdas 9, 40, 5, 4, 119174619Sdas 11, 60, 6, 5, 120174619Sdas 13, 84, 7, 6, 121174619Sdas 33, 56, 7, 4, 122174619Sdas 39, 80, 8, 5, 123174619Sdas 65, 72, 9, 4, 124174619Sdas 987, 9916, 74, 67, 125174619Sdas 5289, 6640, 83, 40, 126174619Sdas 460766389075.0, 16762287900.0, 678910, 12345 127174619Sdas }; 128174619Sdas /* 129174619Sdas * We also test some multiples of the above arguments. This 130174619Sdas * array defines which multiples we use. Note that these have 131174619Sdas * to be small enough to not cause overflow for float precision 132174619Sdas * with all of the constants in the above table. 133174619Sdas */ 134174619Sdas static const double mults[] = { 135174619Sdas 1, 136174619Sdas 2, 137174619Sdas 3, 138174619Sdas 13, 139174619Sdas 16, 140174619Sdas 0x1.p30, 141174619Sdas 0x1.p-30, 142174619Sdas }; 143174619Sdas 144174619Sdas double a, b; 145174619Sdas double x, y; 146174619Sdas int i, j; 147174619Sdas 148174619Sdas for (i = 0; i < N(tests); i += 4) { 149174619Sdas for (j = 0; j < N(mults); j++) { 150174619Sdas a = tests[i] * mults[j] * mults[j]; 151174619Sdas b = tests[i + 1] * mults[j] * mults[j]; 152174619Sdas x = tests[i + 2] * mults[j]; 153174619Sdas y = tests[i + 3] * mults[j]; 154174619Sdas assert(t_csqrt(cpack(a, b)) == cpack(x, y)); 155174619Sdas } 156174619Sdas } 157174619Sdas 158174619Sdas} 159174619Sdas 160174619Sdas/* 161174619Sdas * Test the handling of +/- 0. 162174619Sdas */ 163174619Sdasstatic void 164174619Sdastest_zeros() 165174619Sdas{ 166174619Sdas 167174619Sdas assert_equal(t_csqrt(cpack(0.0, 0.0)), cpack(0.0, 0.0)); 168174619Sdas assert_equal(t_csqrt(cpack(-0.0, 0.0)), cpack(0.0, 0.0)); 169174619Sdas assert_equal(t_csqrt(cpack(0.0, -0.0)), cpack(0.0, -0.0)); 170174619Sdas assert_equal(t_csqrt(cpack(-0.0, -0.0)), cpack(0.0, -0.0)); 171174619Sdas} 172174619Sdas 173174619Sdas/* 174174619Sdas * Test the handling of infinities when the other argument is not NaN. 175174619Sdas */ 176174619Sdasstatic void 177174619Sdastest_infinities() 178174619Sdas{ 179174619Sdas static const double vals[] = { 180174619Sdas 0.0, 181174619Sdas -0.0, 182174619Sdas 42.0, 183174619Sdas -42.0, 184174619Sdas INFINITY, 185174619Sdas -INFINITY, 186174619Sdas }; 187174619Sdas 188174619Sdas int i; 189174619Sdas 190174619Sdas for (i = 0; i < N(vals); i++) { 191174619Sdas if (isfinite(vals[i])) { 192174619Sdas assert_equal(t_csqrt(cpack(-INFINITY, vals[i])), 193174619Sdas cpack(0.0, copysign(INFINITY, vals[i]))); 194174619Sdas assert_equal(t_csqrt(cpack(INFINITY, vals[i])), 195174619Sdas cpack(INFINITY, copysign(0.0, vals[i]))); 196174619Sdas } 197174619Sdas assert_equal(t_csqrt(cpack(vals[i], INFINITY)), 198174619Sdas cpack(INFINITY, INFINITY)); 199174619Sdas assert_equal(t_csqrt(cpack(vals[i], -INFINITY)), 200174619Sdas cpack(INFINITY, -INFINITY)); 201174619Sdas } 202174619Sdas} 203174619Sdas 204174619Sdas/* 205174619Sdas * Test the handling of NaNs. 206174619Sdas */ 207174619Sdasstatic void 208174619Sdastest_nans() 209174619Sdas{ 210174619Sdas 211174619Sdas assert(creal(t_csqrt(cpack(INFINITY, NAN))) == INFINITY); 212174619Sdas assert(isnan(cimag(t_csqrt(cpack(INFINITY, NAN))))); 213174619Sdas 214174619Sdas assert(isnan(creal(t_csqrt(cpack(-INFINITY, NAN))))); 215174619Sdas assert(isinf(cimag(t_csqrt(cpack(-INFINITY, NAN))))); 216174619Sdas 217174619Sdas assert_equal(t_csqrt(cpack(NAN, INFINITY)), cpack(INFINITY, INFINITY)); 218174619Sdas assert_equal(t_csqrt(cpack(NAN, -INFINITY)), 219174619Sdas cpack(INFINITY, -INFINITY)); 220174619Sdas 221174619Sdas assert_equal(t_csqrt(cpack(0.0, NAN)), cpack(NAN, NAN)); 222174619Sdas assert_equal(t_csqrt(cpack(-0.0, NAN)), cpack(NAN, NAN)); 223174619Sdas assert_equal(t_csqrt(cpack(42.0, NAN)), cpack(NAN, NAN)); 224174619Sdas assert_equal(t_csqrt(cpack(-42.0, NAN)), cpack(NAN, NAN)); 225174619Sdas assert_equal(t_csqrt(cpack(NAN, 0.0)), cpack(NAN, NAN)); 226174619Sdas assert_equal(t_csqrt(cpack(NAN, -0.0)), cpack(NAN, NAN)); 227174619Sdas assert_equal(t_csqrt(cpack(NAN, 42.0)), cpack(NAN, NAN)); 228174619Sdas assert_equal(t_csqrt(cpack(NAN, -42.0)), cpack(NAN, NAN)); 229174619Sdas assert_equal(t_csqrt(cpack(NAN, NAN)), cpack(NAN, NAN)); 230174619Sdas} 231174619Sdas 232174619Sdas/* 233174619Sdas * Test whether csqrt(a + bi) works for inputs that are large enough to 234174619Sdas * cause overflow in hypot(a, b) + a. In this case we are using 235174619Sdas * csqrt(115 + 252*I) == 14 + 9*I 236174619Sdas * scaled up to near MAX_EXP. 237174619Sdas */ 238174619Sdasstatic void 239174619Sdastest_overflow(int maxexp) 240174619Sdas{ 241174619Sdas double a, b; 242174619Sdas double complex result; 243174619Sdas 244174619Sdas a = ldexp(115 * 0x1p-8, maxexp); 245174619Sdas b = ldexp(252 * 0x1p-8, maxexp); 246174619Sdas result = t_csqrt(cpack(a, b)); 247174619Sdas assert(creal(result) == ldexp(14 * 0x1p-4, maxexp / 2)); 248174619Sdas assert(cimag(result) == ldexp(9 * 0x1p-4, maxexp / 2)); 249174619Sdas} 250174619Sdas 251174619Sdasint 252174619Sdasmain(int argc, char *argv[]) 253174619Sdas{ 254174619Sdas 255174619Sdas printf("1..10\n"); 256174619Sdas 257174619Sdas /* Test csqrt() */ 258174619Sdas t_csqrt = csqrt; 259174619Sdas 260174619Sdas test_finite(); 261174619Sdas printf("ok 1 - csqrt\n"); 262174619Sdas 263174619Sdas test_zeros(); 264174619Sdas printf("ok 2 - csqrt\n"); 265174619Sdas 266174619Sdas test_infinities(); 267174619Sdas printf("ok 3 - csqrt\n"); 268174619Sdas 269174619Sdas test_nans(); 270174619Sdas printf("ok 4 - csqrt\n"); 271174619Sdas 272174619Sdas test_overflow(DBL_MAX_EXP); 273174619Sdas printf("ok 5 - csqrt\n"); 274174619Sdas 275174619Sdas /* Now test csqrtf() */ 276174619Sdas t_csqrt = _csqrtf; 277174619Sdas 278174619Sdas test_finite(); 279174619Sdas printf("ok 6 - csqrt\n"); 280174619Sdas 281174619Sdas test_zeros(); 282174619Sdas printf("ok 7 - csqrt\n"); 283174619Sdas 284174619Sdas test_infinities(); 285174619Sdas printf("ok 8 - csqrt\n"); 286174619Sdas 287174619Sdas test_nans(); 288174619Sdas printf("ok 9 - csqrt\n"); 289174619Sdas 290174619Sdas test_overflow(FLT_MAX_EXP); 291174619Sdas printf("ok 10 - csqrt\n"); 292174619Sdas 293174619Sdas return (0); 294174619Sdas} 295