1 /* $OpenBSD: cexp_test.c,v 1.1 2021/10/22 18:00:22 mbuhl Exp $ */ 2 /*- 3 * Copyright (c) 2008-2011 David Schultz <das@FreeBSD.org> 4 * All rights reserved. 5 * 6 * Redistribution and use in source and binary forms, with or without 7 * modification, are permitted provided that the following conditions 8 * are met: 9 * 1. Redistributions of source code must retain the above copyright 10 * notice, this list of conditions and the following disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 18 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 19 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 20 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 21 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 25 * SUCH DAMAGE. 26 */ 27 28 #include "macros.h" 29 30 /* 31 * Tests for corner cases in cexp*(). 32 */ 33 34 #include <sys/cdefs.h> 35 __FBSDID("$FreeBSD$"); 36 37 #include <sys/param.h> 38 39 #include <complex.h> 40 #include <fenv.h> 41 #include <float.h> 42 #include <math.h> 43 #include <stdio.h> 44 45 #include "test-utils.h" 46 47 #pragma STDC FENV_ACCESS ON 48 #pragma STDC CX_LIMITED_RANGE OFF 49 50 /* 51 * Test that a function returns the correct value and sets the 52 * exception flags correctly. The exceptmask specifies which 53 * exceptions we should check. We need to be lenient for several 54 * reasons, but mainly because on some architectures it's impossible 55 * to raise FE_OVERFLOW without raising FE_INEXACT. In some cases, 56 * whether cexp() raises an invalid exception is unspecified. 57 * 58 * These are macros instead of functions so that assert provides more 59 * meaningful error messages. 60 * 61 * XXX The volatile here is to avoid gcc's bogus constant folding and work 62 * around the lack of support for the FENV_ACCESS pragma. 63 */ 64 #define test_t(type, func, z, result, exceptmask, excepts, checksign) \ 65 do { \ 66 volatile long double complex _d = z; \ 67 volatile type complex _r = result; \ 68 ATF_REQUIRE_EQ(0, feclearexcept(FE_ALL_EXCEPT)); \ 69 CHECK_CFPEQUAL_CS((func)(_d), (_r), (checksign)); \ 70 CHECK_FP_EXCEPTIONS_MSG(excepts, exceptmask, "for %s(%s)", \ 71 #func, #z); \ 72 } while (0) 73 74 #define test(func, z, result, exceptmask, excepts, checksign) \ 75 test_t(double, func, z, result, exceptmask, excepts, checksign) 76 77 #define test_f(func, z, result, exceptmask, excepts, checksign) \ 78 test_t(float, func, z, result, exceptmask, excepts, checksign) 79 80 #define test_l(func, z, result, exceptmask, excepts, checksign) \ 81 test_t(long double, func, z, result, exceptmask, excepts, \ 82 checksign) 83 /* Test within a given tolerance. */ 84 #define test_tol(func, z, result, tol) do { \ 85 CHECK_CFPEQUAL_TOL((func)(z), (result), (tol), \ 86 FPE_ABS_ZERO | CS_BOTH); \ 87 } while (0) 88 89 /* Test all the functions that compute cexp(x). */ 90 #define testall(x, result, exceptmask, excepts, checksign) do { \ 91 test(cexp, x, result, exceptmask, excepts, checksign); \ 92 test_f(cexpf, x, result, exceptmask, excepts, checksign); \ 93 test_l(cexpl, x, result, exceptmask, excepts, checksign); \ 94 } while (0) 95 96 /* 97 * Test all the functions that compute cexp(x), within a given tolerance. 98 * The tolerance is specified in ulps. 99 */ 100 #define testall_tol(x, result, tol) do { \ 101 test_tol(cexp, x, result, tol * DBL_ULP()); \ 102 test_tol(cexpf, x, result, tol * FLT_ULP()); \ 103 } while (0) 104 105 /* Various finite non-zero numbers to test. */ 106 static const float finites[] = 107 { -42.0e20, -1.0, -1.0e-10, -0.0, 0.0, 1.0e-10, 1.0, 42.0e20 }; 108 109 110 /* Tests for 0 */ 111 ATF_TC_WITHOUT_HEAD(zero); 112 ATF_TC_BODY(zero, tc) 113 { 114 115 /* cexp(0) = 1, no exceptions raised */ 116 testall(0.0, 1.0, ALL_STD_EXCEPT, 0, 1); 117 testall(-0.0, 1.0, ALL_STD_EXCEPT, 0, 1); 118 testall(CMPLXL(0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1); 119 testall(CMPLXL(-0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1); 120 } 121 122 /* 123 * Tests for NaN. The signs of the results are indeterminate unless the 124 * imaginary part is 0. 125 */ 126 ATF_TC_WITHOUT_HEAD(nan); 127 ATF_TC_BODY(nan, tc) 128 { 129 unsigned i; 130 131 /* cexp(x + NaNi) = NaN + NaNi and optionally raises invalid */ 132 /* cexp(NaN + yi) = NaN + NaNi and optionally raises invalid (|y|>0) */ 133 for (i = 0; i < nitems(finites); i++) { 134 testall(CMPLXL(finites[i], NAN), CMPLXL(NAN, NAN), 135 ALL_STD_EXCEPT & ~FE_INVALID, 0, 0); 136 if (finites[i] == 0.0) 137 continue; 138 #ifndef __OpenBSD__ 139 /* XXX FE_INEXACT shouldn't be raised here */ 140 testall(CMPLXL(NAN, finites[i]), CMPLXL(NAN, NAN), 141 ALL_STD_EXCEPT & ~(FE_INVALID | FE_INEXACT), 0, 0); 142 #else 143 testall(CMPLXL(NAN, finites[i]), CMPLXL(NAN, NAN), 144 ALL_STD_EXCEPT & ~(FE_INVALID), 0, 0); 145 #endif 146 } 147 148 /* cexp(NaN +- 0i) = NaN +- 0i */ 149 testall(CMPLXL(NAN, 0.0), CMPLXL(NAN, 0.0), ALL_STD_EXCEPT, 0, 1); 150 testall(CMPLXL(NAN, -0.0), CMPLXL(NAN, -0.0), ALL_STD_EXCEPT, 0, 1); 151 152 /* cexp(inf + NaN i) = inf + nan i */ 153 testall(CMPLXL(INFINITY, NAN), CMPLXL(INFINITY, NAN), 154 ALL_STD_EXCEPT, 0, 0); 155 /* cexp(-inf + NaN i) = 0 */ 156 testall(CMPLXL(-INFINITY, NAN), CMPLXL(0.0, 0.0), 157 ALL_STD_EXCEPT, 0, 0); 158 /* cexp(NaN + NaN i) = NaN + NaN i */ 159 testall(CMPLXL(NAN, NAN), CMPLXL(NAN, NAN), 160 ALL_STD_EXCEPT, 0, 0); 161 } 162 163 ATF_TC_WITHOUT_HEAD(inf); 164 ATF_TC_BODY(inf, tc) 165 { 166 unsigned i; 167 168 /* cexp(x + inf i) = NaN + NaNi and raises invalid */ 169 for (i = 0; i < nitems(finites); i++) { 170 testall(CMPLXL(finites[i], INFINITY), CMPLXL(NAN, NAN), 171 ALL_STD_EXCEPT, FE_INVALID, 1); 172 } 173 /* cexp(-inf + yi) = 0 * (cos(y) + sin(y)i) */ 174 /* XXX shouldn't raise an inexact exception */ 175 testall(CMPLXL(-INFINITY, M_PI_4), CMPLXL(0.0, 0.0), 176 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 177 testall(CMPLXL(-INFINITY, 3 * M_PI_4), CMPLXL(-0.0, 0.0), 178 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 179 testall(CMPLXL(-INFINITY, 5 * M_PI_4), CMPLXL(-0.0, -0.0), 180 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 181 testall(CMPLXL(-INFINITY, 7 * M_PI_4), CMPLXL(0.0, -0.0), 182 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 183 testall(CMPLXL(-INFINITY, 0.0), CMPLXL(0.0, 0.0), 184 ALL_STD_EXCEPT, 0, 1); 185 testall(CMPLXL(-INFINITY, -0.0), CMPLXL(0.0, -0.0), 186 ALL_STD_EXCEPT, 0, 1); 187 /* cexp(inf + yi) = inf * (cos(y) + sin(y)i) (except y=0) */ 188 /* XXX shouldn't raise an inexact exception */ 189 testall(CMPLXL(INFINITY, M_PI_4), CMPLXL(INFINITY, INFINITY), 190 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 191 testall(CMPLXL(INFINITY, 3 * M_PI_4), CMPLXL(-INFINITY, INFINITY), 192 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 193 testall(CMPLXL(INFINITY, 5 * M_PI_4), CMPLXL(-INFINITY, -INFINITY), 194 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 195 testall(CMPLXL(INFINITY, 7 * M_PI_4), CMPLXL(INFINITY, -INFINITY), 196 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 197 /* cexp(inf + 0i) = inf + 0i */ 198 testall(CMPLXL(INFINITY, 0.0), CMPLXL(INFINITY, 0.0), 199 ALL_STD_EXCEPT, 0, 1); 200 testall(CMPLXL(INFINITY, -0.0), CMPLXL(INFINITY, -0.0), 201 ALL_STD_EXCEPT, 0, 1); 202 } 203 204 ATF_TC_WITHOUT_HEAD(reals); 205 ATF_TC_BODY(reals, tc) 206 { 207 unsigned i; 208 209 for (i = 0; i < nitems(finites); i++) { 210 /* XXX could check exceptions more meticulously */ 211 test(cexp, CMPLXL(finites[i], 0.0), 212 CMPLXL(exp(finites[i]), 0.0), 213 FE_INVALID | FE_DIVBYZERO, 0, 1); 214 test(cexp, CMPLXL(finites[i], -0.0), 215 CMPLXL(exp(finites[i]), -0.0), 216 FE_INVALID | FE_DIVBYZERO, 0, 1); 217 test_f(cexpf, CMPLXL(finites[i], 0.0), 218 CMPLXL(expf(finites[i]), 0.0), 219 FE_INVALID | FE_DIVBYZERO, 0, 1); 220 test_f(cexpf, CMPLXL(finites[i], -0.0), 221 CMPLXL(expf(finites[i]), -0.0), 222 FE_INVALID | FE_DIVBYZERO, 0, 1); 223 } 224 } 225 226 ATF_TC_WITHOUT_HEAD(imaginaries); 227 ATF_TC_BODY(imaginaries, tc) 228 { 229 unsigned i; 230 231 for (i = 0; i < nitems(finites); i++) { 232 test(cexp, CMPLXL(0.0, finites[i]), 233 CMPLXL(cos(finites[i]), sin(finites[i])), 234 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 235 test(cexp, CMPLXL(-0.0, finites[i]), 236 CMPLXL(cos(finites[i]), sin(finites[i])), 237 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 238 test_f(cexpf, CMPLXL(0.0, finites[i]), 239 CMPLXL(cosf(finites[i]), sinf(finites[i])), 240 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 241 test_f(cexpf, CMPLXL(-0.0, finites[i]), 242 CMPLXL(cosf(finites[i]), sinf(finites[i])), 243 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 244 } 245 } 246 247 ATF_TC_WITHOUT_HEAD(small); 248 ATF_TC_BODY(small, tc) 249 { 250 static const double tests[] = { 251 /* csqrt(a + bI) = x + yI */ 252 /* a b x y */ 253 1.0, M_PI_4, M_SQRT2 * 0.5 * M_E, M_SQRT2 * 0.5 * M_E, 254 -1.0, M_PI_4, M_SQRT2 * 0.5 / M_E, M_SQRT2 * 0.5 / M_E, 255 2.0, M_PI_2, 0.0, M_E * M_E, 256 M_LN2, M_PI, -2.0, 0.0, 257 }; 258 double a, b; 259 double x, y; 260 unsigned i; 261 262 for (i = 0; i < nitems(tests); i += 4) { 263 a = tests[i]; 264 b = tests[i + 1]; 265 x = tests[i + 2]; 266 y = tests[i + 3]; 267 test_tol(cexp, CMPLXL(a, b), CMPLXL(x, y), 3 * DBL_ULP()); 268 269 /* float doesn't have enough precision to pass these tests */ 270 if (x == 0 || y == 0) 271 continue; 272 test_tol(cexpf, CMPLXL(a, b), CMPLXL(x, y), 1 * FLT_ULP()); 273 } 274 } 275 276 /* Test inputs with a real part r that would overflow exp(r). */ 277 ATF_TC_WITHOUT_HEAD(large); 278 ATF_TC_BODY(large, tc) 279 { 280 281 test_tol(cexp, CMPLXL(709.79, 0x1p-1074), 282 CMPLXL(INFINITY, 8.94674309915433533273e-16), DBL_ULP()); 283 test_tol(cexp, CMPLXL(1000, 0x1p-1074), 284 CMPLXL(INFINITY, 9.73344457300016401328e+110), DBL_ULP()); 285 test_tol(cexp, CMPLXL(1400, 0x1p-1074), 286 CMPLXL(INFINITY, 5.08228858149196559681e+284), DBL_ULP()); 287 test_tol(cexp, CMPLXL(900, 0x1.23456789abcdep-1020), 288 CMPLXL(INFINITY, 7.42156649354218408074e+83), DBL_ULP()); 289 test_tol(cexp, CMPLXL(1300, 0x1.23456789abcdep-1020), 290 CMPLXL(INFINITY, 3.87514844965996756704e+257), DBL_ULP()); 291 292 test_tol(cexpf, CMPLXL(88.73, 0x1p-149), 293 CMPLXL(INFINITY, 4.80265603e-07), 2 * FLT_ULP()); 294 test_tol(cexpf, CMPLXL(90, 0x1p-149), 295 CMPLXL(INFINITY, 1.7101492622e-06f), 2 * FLT_ULP()); 296 test_tol(cexpf, CMPLXL(192, 0x1p-149), 297 CMPLXL(INFINITY, 3.396809344e+38f), 2 * FLT_ULP()); 298 test_tol(cexpf, CMPLXL(120, 0x1.234568p-120), 299 CMPLXL(INFINITY, 1.1163382522e+16f), 2 * FLT_ULP()); 300 test_tol(cexpf, CMPLXL(170, 0x1.234568p-120), 301 CMPLXL(INFINITY, 5.7878851079e+37f), 2 * FLT_ULP()); 302 } 303 304 ATF_TP_ADD_TCS(tp) 305 { 306 ATF_TP_ADD_TC(tp, zero); 307 ATF_TP_ADD_TC(tp, nan); 308 ATF_TP_ADD_TC(tp, inf); 309 ATF_TP_ADD_TC(tp, reals); 310 ATF_TP_ADD_TC(tp, imaginaries); 311 ATF_TP_ADD_TC(tp, small); 312 ATF_TP_ADD_TC(tp, large); 313 314 return (atf_no_error()); 315 } 316