csqrt_test.c revision 177763
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 177763 2008-03-30 20:09:51Z 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/* 43177763Sdas * This is a test hook that can point to csqrtl(), _csqrt(), or to _csqrtf(). 44177763Sdas * The latter two convert to float or double, respectively, and test csqrtf() 45177763Sdas * and csqrt() with the same arguments. 46174619Sdas */ 47177763Sdaslong double complex (*t_csqrt)(long double complex); 48174619Sdas 49177763Sdasstatic long double complex 50177763Sdas_csqrtf(long double complex d) 51174619Sdas{ 52174619Sdas 53174619Sdas return (csqrtf((float complex)d)); 54174619Sdas} 55174619Sdas 56177763Sdasstatic long double complex 57177763Sdas_csqrt(long double complex d) 58177763Sdas{ 59177763Sdas 60177763Sdas return (csqrt((double complex)d)); 61177763Sdas} 62177763Sdas 63174619Sdas#pragma STDC CX_LIMITED_RANGE off 64174619Sdas 65174619Sdas/* 66174619Sdas * XXX gcc implements complex multiplication incorrectly. In 67174619Sdas * particular, it implements it as if the CX_LIMITED_RANGE pragma 68174619Sdas * were ON. Consequently, we need this function to form numbers 69174619Sdas * such as x + INFINITY * I, since gcc evalutes INFINITY * I as 70174619Sdas * NaN + INFINITY * I. 71174619Sdas */ 72177763Sdasstatic inline long double complex 73177763Sdascpackl(long double x, long double y) 74174619Sdas{ 75177763Sdas long double complex z; 76174619Sdas 77174619Sdas __real__ z = x; 78174619Sdas __imag__ z = y; 79174619Sdas return (z); 80174619Sdas} 81174619Sdas 82174619Sdas/* 83174619Sdas * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0. 84174619Sdas * Fail an assertion if they differ. 85174619Sdas */ 86174619Sdasstatic void 87177763Sdasassert_equal(long double complex d1, long double complex d2) 88174619Sdas{ 89174619Sdas 90177763Sdas if (isnan(creall(d1))) { 91177763Sdas assert(isnan(creall(d2))); 92174619Sdas } else { 93177763Sdas assert(creall(d1) == creall(d2)); 94177763Sdas assert(copysignl(1.0, creall(d1)) == 95177763Sdas copysignl(1.0, creall(d2))); 96174619Sdas } 97177763Sdas if (isnan(cimagl(d1))) { 98177763Sdas assert(isnan(cimagl(d2))); 99174619Sdas } else { 100177763Sdas assert(cimagl(d1) == cimagl(d2)); 101177763Sdas assert(copysignl(1.0, cimagl(d1)) == 102177763Sdas copysignl(1.0, cimagl(d2))); 103174619Sdas } 104174619Sdas} 105174619Sdas 106174619Sdas/* 107174619Sdas * Test csqrt for some finite arguments where the answer is exact. 108174619Sdas * (We do not test if it produces correctly rounded answers when the 109174619Sdas * result is inexact, nor do we check whether it throws spurious 110174619Sdas * exceptions.) 111174619Sdas */ 112174619Sdasstatic void 113174619Sdastest_finite() 114174619Sdas{ 115174619Sdas static const double tests[] = { 116174619Sdas /* csqrt(a + bI) = x + yI */ 117174619Sdas /* a b x y */ 118174619Sdas 0, 8, 2, 2, 119174619Sdas 0, -8, 2, -2, 120174619Sdas 4, 0, 2, 0, 121174619Sdas -4, 0, 0, 2, 122174619Sdas 3, 4, 2, 1, 123174619Sdas 3, -4, 2, -1, 124174619Sdas -3, 4, 1, 2, 125174619Sdas -3, -4, 1, -2, 126174619Sdas 5, 12, 3, 2, 127174619Sdas 7, 24, 4, 3, 128174619Sdas 9, 40, 5, 4, 129174619Sdas 11, 60, 6, 5, 130174619Sdas 13, 84, 7, 6, 131174619Sdas 33, 56, 7, 4, 132174619Sdas 39, 80, 8, 5, 133174619Sdas 65, 72, 9, 4, 134174619Sdas 987, 9916, 74, 67, 135174619Sdas 5289, 6640, 83, 40, 136174619Sdas 460766389075.0, 16762287900.0, 678910, 12345 137174619Sdas }; 138174619Sdas /* 139174619Sdas * We also test some multiples of the above arguments. This 140174619Sdas * array defines which multiples we use. Note that these have 141174619Sdas * to be small enough to not cause overflow for float precision 142174619Sdas * with all of the constants in the above table. 143174619Sdas */ 144174619Sdas static const double mults[] = { 145174619Sdas 1, 146174619Sdas 2, 147174619Sdas 3, 148174619Sdas 13, 149174619Sdas 16, 150174619Sdas 0x1.p30, 151174619Sdas 0x1.p-30, 152174619Sdas }; 153174619Sdas 154174619Sdas double a, b; 155174619Sdas double x, y; 156174619Sdas int i, j; 157174619Sdas 158174619Sdas for (i = 0; i < N(tests); i += 4) { 159174619Sdas for (j = 0; j < N(mults); j++) { 160174619Sdas a = tests[i] * mults[j] * mults[j]; 161174619Sdas b = tests[i + 1] * mults[j] * mults[j]; 162174619Sdas x = tests[i + 2] * mults[j]; 163174619Sdas y = tests[i + 3] * mults[j]; 164177763Sdas assert(t_csqrt(cpackl(a, b)) == cpackl(x, y)); 165174619Sdas } 166174619Sdas } 167174619Sdas 168174619Sdas} 169174619Sdas 170174619Sdas/* 171174619Sdas * Test the handling of +/- 0. 172174619Sdas */ 173174619Sdasstatic void 174174619Sdastest_zeros() 175174619Sdas{ 176174619Sdas 177177763Sdas assert_equal(t_csqrt(cpackl(0.0, 0.0)), cpackl(0.0, 0.0)); 178177763Sdas assert_equal(t_csqrt(cpackl(-0.0, 0.0)), cpackl(0.0, 0.0)); 179177763Sdas assert_equal(t_csqrt(cpackl(0.0, -0.0)), cpackl(0.0, -0.0)); 180177763Sdas assert_equal(t_csqrt(cpackl(-0.0, -0.0)), cpackl(0.0, -0.0)); 181174619Sdas} 182174619Sdas 183174619Sdas/* 184174619Sdas * Test the handling of infinities when the other argument is not NaN. 185174619Sdas */ 186174619Sdasstatic void 187174619Sdastest_infinities() 188174619Sdas{ 189174619Sdas static const double vals[] = { 190174619Sdas 0.0, 191174619Sdas -0.0, 192174619Sdas 42.0, 193174619Sdas -42.0, 194174619Sdas INFINITY, 195174619Sdas -INFINITY, 196174619Sdas }; 197174619Sdas 198174619Sdas int i; 199174619Sdas 200174619Sdas for (i = 0; i < N(vals); i++) { 201174619Sdas if (isfinite(vals[i])) { 202177763Sdas assert_equal(t_csqrt(cpackl(-INFINITY, vals[i])), 203177763Sdas cpackl(0.0, copysignl(INFINITY, vals[i]))); 204177763Sdas assert_equal(t_csqrt(cpackl(INFINITY, vals[i])), 205177763Sdas cpackl(INFINITY, copysignl(0.0, vals[i]))); 206174619Sdas } 207177763Sdas assert_equal(t_csqrt(cpackl(vals[i], INFINITY)), 208177763Sdas cpackl(INFINITY, INFINITY)); 209177763Sdas assert_equal(t_csqrt(cpackl(vals[i], -INFINITY)), 210177763Sdas cpackl(INFINITY, -INFINITY)); 211174619Sdas } 212174619Sdas} 213174619Sdas 214174619Sdas/* 215174619Sdas * Test the handling of NaNs. 216174619Sdas */ 217174619Sdasstatic void 218174619Sdastest_nans() 219174619Sdas{ 220174619Sdas 221177763Sdas assert(creall(t_csqrt(cpackl(INFINITY, NAN))) == INFINITY); 222177763Sdas assert(isnan(cimagl(t_csqrt(cpackl(INFINITY, NAN))))); 223174619Sdas 224177763Sdas assert(isnan(creall(t_csqrt(cpackl(-INFINITY, NAN))))); 225177763Sdas assert(isinf(cimagl(t_csqrt(cpackl(-INFINITY, NAN))))); 226174619Sdas 227177763Sdas assert_equal(t_csqrt(cpackl(NAN, INFINITY)), 228177763Sdas cpackl(INFINITY, INFINITY)); 229177763Sdas assert_equal(t_csqrt(cpackl(NAN, -INFINITY)), 230177763Sdas cpackl(INFINITY, -INFINITY)); 231174619Sdas 232177763Sdas assert_equal(t_csqrt(cpackl(0.0, NAN)), cpackl(NAN, NAN)); 233177763Sdas assert_equal(t_csqrt(cpackl(-0.0, NAN)), cpackl(NAN, NAN)); 234177763Sdas assert_equal(t_csqrt(cpackl(42.0, NAN)), cpackl(NAN, NAN)); 235177763Sdas assert_equal(t_csqrt(cpackl(-42.0, NAN)), cpackl(NAN, NAN)); 236177763Sdas assert_equal(t_csqrt(cpackl(NAN, 0.0)), cpackl(NAN, NAN)); 237177763Sdas assert_equal(t_csqrt(cpackl(NAN, -0.0)), cpackl(NAN, NAN)); 238177763Sdas assert_equal(t_csqrt(cpackl(NAN, 42.0)), cpackl(NAN, NAN)); 239177763Sdas assert_equal(t_csqrt(cpackl(NAN, -42.0)), cpackl(NAN, NAN)); 240177763Sdas assert_equal(t_csqrt(cpackl(NAN, NAN)), cpackl(NAN, NAN)); 241174619Sdas} 242174619Sdas 243174619Sdas/* 244174619Sdas * Test whether csqrt(a + bi) works for inputs that are large enough to 245174619Sdas * cause overflow in hypot(a, b) + a. In this case we are using 246174619Sdas * csqrt(115 + 252*I) == 14 + 9*I 247174619Sdas * scaled up to near MAX_EXP. 248174619Sdas */ 249174619Sdasstatic void 250174619Sdastest_overflow(int maxexp) 251174619Sdas{ 252177763Sdas long double a, b; 253177763Sdas long double complex result; 254174619Sdas 255177763Sdas a = ldexpl(115 * 0x1p-8, maxexp); 256177763Sdas b = ldexpl(252 * 0x1p-8, maxexp); 257177763Sdas result = t_csqrt(cpackl(a, b)); 258177763Sdas assert(creall(result) == ldexpl(14 * 0x1p-4, maxexp / 2)); 259177763Sdas assert(cimagl(result) == ldexpl(9 * 0x1p-4, maxexp / 2)); 260174619Sdas} 261174619Sdas 262174619Sdasint 263174619Sdasmain(int argc, char *argv[]) 264174619Sdas{ 265174619Sdas 266177763Sdas printf("1..15\n"); 267174619Sdas 268174619Sdas /* Test csqrt() */ 269177763Sdas t_csqrt = _csqrt; 270174619Sdas 271174619Sdas test_finite(); 272174619Sdas printf("ok 1 - csqrt\n"); 273174619Sdas 274174619Sdas test_zeros(); 275174619Sdas printf("ok 2 - csqrt\n"); 276174619Sdas 277174619Sdas test_infinities(); 278174619Sdas printf("ok 3 - csqrt\n"); 279174619Sdas 280174619Sdas test_nans(); 281174619Sdas printf("ok 4 - csqrt\n"); 282174619Sdas 283174619Sdas test_overflow(DBL_MAX_EXP); 284174619Sdas printf("ok 5 - csqrt\n"); 285174619Sdas 286174619Sdas /* Now test csqrtf() */ 287174619Sdas t_csqrt = _csqrtf; 288174619Sdas 289174619Sdas test_finite(); 290174619Sdas printf("ok 6 - csqrt\n"); 291174619Sdas 292174619Sdas test_zeros(); 293174619Sdas printf("ok 7 - csqrt\n"); 294174619Sdas 295174619Sdas test_infinities(); 296174619Sdas printf("ok 8 - csqrt\n"); 297174619Sdas 298174619Sdas test_nans(); 299174619Sdas printf("ok 9 - csqrt\n"); 300174619Sdas 301174619Sdas test_overflow(FLT_MAX_EXP); 302174619Sdas printf("ok 10 - csqrt\n"); 303174619Sdas 304177763Sdas /* Now test csqrtl() */ 305177763Sdas t_csqrt = csqrtl; 306177763Sdas 307177763Sdas test_finite(); 308177763Sdas printf("ok 11 - csqrt\n"); 309177763Sdas 310177763Sdas test_zeros(); 311177763Sdas printf("ok 12 - csqrt\n"); 312177763Sdas 313177763Sdas test_infinities(); 314177763Sdas printf("ok 13 - csqrt\n"); 315177763Sdas 316177763Sdas test_nans(); 317177763Sdas printf("ok 14 - csqrt\n"); 318177763Sdas 319177763Sdas test_overflow(LDBL_MAX_EXP); 320177763Sdas printf("ok 15 - csqrt\n"); 321177763Sdas 322174619Sdas return (0); 323174619Sdas} 324