1/* $OpenBSD: cexp_test.c,v 1.3 2021/12/13 18:04:28 deraadt 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/types.h> 35 36#include <complex.h> 37#include <fenv.h> 38#include <float.h> 39#include <math.h> 40#include <stdio.h> 41 42#include "test-utils.h" 43 44#pragma STDC FENV_ACCESS ON 45#pragma STDC CX_LIMITED_RANGE OFF 46 47/* 48 * Test that a function returns the correct value and sets the 49 * exception flags correctly. The exceptmask specifies which 50 * exceptions we should check. We need to be lenient for several 51 * reasons, but mainly because on some architectures it's impossible 52 * to raise FE_OVERFLOW without raising FE_INEXACT. In some cases, 53 * whether cexp() raises an invalid exception is unspecified. 54 * 55 * These are macros instead of functions so that assert provides more 56 * meaningful error messages. 57 * 58 * XXX The volatile here is to avoid gcc's bogus constant folding and work 59 * around the lack of support for the FENV_ACCESS pragma. 60 */ 61#define test_t(type, func, z, result, exceptmask, excepts, checksign) \ 62do { \ 63 volatile long double complex _d = z; \ 64 volatile type complex _r = result; \ 65 ATF_REQUIRE_EQ(0, feclearexcept(FE_ALL_EXCEPT)); \ 66 CHECK_CFPEQUAL_CS((func)(_d), (_r), (checksign)); \ 67 CHECK_FP_EXCEPTIONS_MSG(excepts, exceptmask, "for %s(%s)", \ 68 #func, #z); \ 69} while (0) 70 71#define test(func, z, result, exceptmask, excepts, checksign) \ 72 test_t(double, func, z, result, exceptmask, excepts, checksign) 73 74#define test_f(func, z, result, exceptmask, excepts, checksign) \ 75 test_t(float, func, z, result, exceptmask, excepts, checksign) 76 77#define test_l(func, z, result, exceptmask, excepts, checksign) \ 78 test_t(long double, func, z, result, exceptmask, excepts, \ 79 checksign) 80/* Test within a given tolerance. */ 81#define test_tol(func, z, result, tol) do { \ 82 CHECK_CFPEQUAL_TOL((func)(z), (result), (tol), \ 83 FPE_ABS_ZERO | CS_BOTH); \ 84} while (0) 85 86/* Test all the functions that compute cexp(x). */ 87#define testall(x, result, exceptmask, excepts, checksign) do { \ 88 test(cexp, x, result, exceptmask, excepts, checksign); \ 89 test_f(cexpf, x, result, exceptmask, excepts, checksign); \ 90 test_l(cexpl, x, result, exceptmask, excepts, checksign); \ 91} while (0) 92 93/* 94 * Test all the functions that compute cexp(x), within a given tolerance. 95 * The tolerance is specified in ulps. 96 */ 97#define testall_tol(x, result, tol) do { \ 98 test_tol(cexp, x, result, tol * DBL_ULP()); \ 99 test_tol(cexpf, x, result, tol * FLT_ULP()); \ 100} while (0) 101 102/* Various finite non-zero numbers to test. */ 103static const float finites[] = 104{ -42.0e20, -1.0, -1.0e-10, -0.0, 0.0, 1.0e-10, 1.0, 42.0e20 }; 105 106 107/* Tests for 0 */ 108ATF_TC_WITHOUT_HEAD(zero); 109ATF_TC_BODY(zero, tc) 110{ 111 112 /* cexp(0) = 1, no exceptions raised */ 113 testall(0.0, 1.0, ALL_STD_EXCEPT, 0, 1); 114 testall(-0.0, 1.0, ALL_STD_EXCEPT, 0, 1); 115 testall(CMPLXL(0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1); 116 testall(CMPLXL(-0.0, -0.0), CMPLXL(1.0, -0.0), ALL_STD_EXCEPT, 0, 1); 117} 118 119/* 120 * Tests for NaN. The signs of the results are indeterminate unless the 121 * imaginary part is 0. 122 */ 123ATF_TC_WITHOUT_HEAD(nan); 124ATF_TC_BODY(nan, tc) 125{ 126 unsigned i; 127 128 /* cexp(x + NaNi) = NaN + NaNi and optionally raises invalid */ 129 /* cexp(NaN + yi) = NaN + NaNi and optionally raises invalid (|y|>0) */ 130 for (i = 0; i < nitems(finites); i++) { 131 testall(CMPLXL(finites[i], NAN), CMPLXL(NAN, NAN), 132 ALL_STD_EXCEPT & ~FE_INVALID, 0, 0); 133 if (finites[i] == 0.0) 134 continue; 135#ifndef __OpenBSD__ 136 /* XXX FE_INEXACT shouldn't be raised here */ 137 testall(CMPLXL(NAN, finites[i]), CMPLXL(NAN, NAN), 138 ALL_STD_EXCEPT & ~(FE_INVALID | FE_INEXACT), 0, 0); 139#else 140 testall(CMPLXL(NAN, finites[i]), CMPLXL(NAN, NAN), 141 ALL_STD_EXCEPT & ~(FE_INVALID), 0, 0); 142#endif 143 } 144 145 /* cexp(NaN +- 0i) = NaN +- 0i */ 146 testall(CMPLXL(NAN, 0.0), CMPLXL(NAN, 0.0), ALL_STD_EXCEPT, 0, 1); 147 testall(CMPLXL(NAN, -0.0), CMPLXL(NAN, -0.0), ALL_STD_EXCEPT, 0, 1); 148 149 /* cexp(inf + NaN i) = inf + nan i */ 150 testall(CMPLXL(INFINITY, NAN), CMPLXL(INFINITY, NAN), 151 ALL_STD_EXCEPT, 0, 0); 152 /* cexp(-inf + NaN i) = 0 */ 153 testall(CMPLXL(-INFINITY, NAN), CMPLXL(0.0, 0.0), 154 ALL_STD_EXCEPT, 0, 0); 155 /* cexp(NaN + NaN i) = NaN + NaN i */ 156 testall(CMPLXL(NAN, NAN), CMPLXL(NAN, NAN), 157 ALL_STD_EXCEPT, 0, 0); 158} 159 160ATF_TC_WITHOUT_HEAD(inf); 161ATF_TC_BODY(inf, tc) 162{ 163 unsigned i; 164 165 /* cexp(x + inf i) = NaN + NaNi and raises invalid */ 166 for (i = 0; i < nitems(finites); i++) { 167 testall(CMPLXL(finites[i], INFINITY), CMPLXL(NAN, NAN), 168 ALL_STD_EXCEPT, FE_INVALID, 1); 169 } 170 /* cexp(-inf + yi) = 0 * (cos(y) + sin(y)i) */ 171 /* XXX shouldn't raise an inexact exception */ 172 testall(CMPLXL(-INFINITY, M_PI_4), CMPLXL(0.0, 0.0), 173 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 174 testall(CMPLXL(-INFINITY, 3 * M_PI_4), CMPLXL(-0.0, 0.0), 175 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 176 testall(CMPLXL(-INFINITY, 5 * M_PI_4), CMPLXL(-0.0, -0.0), 177 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 178 testall(CMPLXL(-INFINITY, 7 * M_PI_4), CMPLXL(0.0, -0.0), 179 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 180 testall(CMPLXL(-INFINITY, 0.0), CMPLXL(0.0, 0.0), 181 ALL_STD_EXCEPT, 0, 1); 182 testall(CMPLXL(-INFINITY, -0.0), CMPLXL(0.0, -0.0), 183 ALL_STD_EXCEPT, 0, 1); 184 /* cexp(inf + yi) = inf * (cos(y) + sin(y)i) (except y=0) */ 185 /* XXX shouldn't raise an inexact exception */ 186 testall(CMPLXL(INFINITY, M_PI_4), CMPLXL(INFINITY, INFINITY), 187 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 188 testall(CMPLXL(INFINITY, 3 * M_PI_4), CMPLXL(-INFINITY, INFINITY), 189 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 190 testall(CMPLXL(INFINITY, 5 * M_PI_4), CMPLXL(-INFINITY, -INFINITY), 191 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 192 testall(CMPLXL(INFINITY, 7 * M_PI_4), CMPLXL(INFINITY, -INFINITY), 193 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 194 /* cexp(inf + 0i) = inf + 0i */ 195 testall(CMPLXL(INFINITY, 0.0), CMPLXL(INFINITY, 0.0), 196 ALL_STD_EXCEPT, 0, 1); 197 testall(CMPLXL(INFINITY, -0.0), CMPLXL(INFINITY, -0.0), 198 ALL_STD_EXCEPT, 0, 1); 199} 200 201ATF_TC_WITHOUT_HEAD(reals); 202ATF_TC_BODY(reals, tc) 203{ 204 unsigned i; 205 206 for (i = 0; i < nitems(finites); i++) { 207 /* XXX could check exceptions more meticulously */ 208 test(cexp, CMPLXL(finites[i], 0.0), 209 CMPLXL(exp(finites[i]), 0.0), 210 FE_INVALID | FE_DIVBYZERO, 0, 1); 211 test(cexp, CMPLXL(finites[i], -0.0), 212 CMPLXL(exp(finites[i]), -0.0), 213 FE_INVALID | FE_DIVBYZERO, 0, 1); 214 test_f(cexpf, CMPLXL(finites[i], 0.0), 215 CMPLXL(expf(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 } 221} 222 223ATF_TC_WITHOUT_HEAD(imaginaries); 224ATF_TC_BODY(imaginaries, tc) 225{ 226 unsigned i; 227 228 for (i = 0; i < nitems(finites); i++) { 229 test(cexp, CMPLXL(0.0, finites[i]), 230 CMPLXL(cos(finites[i]), sin(finites[i])), 231 ALL_STD_EXCEPT & ~FE_INEXACT, 0, 1); 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_f(cexpf, CMPLXL(0.0, finites[i]), 236 CMPLXL(cosf(finites[i]), sinf(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 } 242} 243 244ATF_TC_WITHOUT_HEAD(small); 245ATF_TC_BODY(small, tc) 246{ 247 static const double tests[] = { 248 /* csqrt(a + bI) = x + yI */ 249 /* a b x y */ 250 1.0, M_PI_4, M_SQRT2 * 0.5 * M_E, M_SQRT2 * 0.5 * M_E, 251 -1.0, M_PI_4, M_SQRT2 * 0.5 / M_E, M_SQRT2 * 0.5 / M_E, 252 2.0, M_PI_2, 0.0, M_E * M_E, 253 M_LN2, M_PI, -2.0, 0.0, 254 }; 255 double a, b; 256 double x, y; 257 unsigned i; 258 259 for (i = 0; i < nitems(tests); i += 4) { 260 a = tests[i]; 261 b = tests[i + 1]; 262 x = tests[i + 2]; 263 y = tests[i + 3]; 264 test_tol(cexp, CMPLXL(a, b), CMPLXL(x, y), 3 * DBL_ULP()); 265 266 /* float doesn't have enough precision to pass these tests */ 267 if (x == 0 || y == 0) 268 continue; 269 test_tol(cexpf, CMPLXL(a, b), CMPLXL(x, y), 1 * FLT_ULP()); 270 } 271} 272 273/* Test inputs with a real part r that would overflow exp(r). */ 274ATF_TC_WITHOUT_HEAD(large); 275ATF_TC_BODY(large, tc) 276{ 277 278 test_tol(cexp, CMPLXL(709.79, 0x1p-1074), 279 CMPLXL(INFINITY, 8.94674309915433533273e-16), DBL_ULP()); 280 test_tol(cexp, CMPLXL(1000, 0x1p-1074), 281 CMPLXL(INFINITY, 9.73344457300016401328e+110), DBL_ULP()); 282 test_tol(cexp, CMPLXL(1400, 0x1p-1074), 283 CMPLXL(INFINITY, 5.08228858149196559681e+284), DBL_ULP()); 284 test_tol(cexp, CMPLXL(900, 0x1.23456789abcdep-1020), 285 CMPLXL(INFINITY, 7.42156649354218408074e+83), DBL_ULP()); 286 test_tol(cexp, CMPLXL(1300, 0x1.23456789abcdep-1020), 287 CMPLXL(INFINITY, 3.87514844965996756704e+257), DBL_ULP()); 288 289 test_tol(cexpf, CMPLXL(88.73, 0x1p-149), 290 CMPLXL(INFINITY, 4.80265603e-07), 2 * FLT_ULP()); 291 test_tol(cexpf, CMPLXL(90, 0x1p-149), 292 CMPLXL(INFINITY, 1.7101492622e-06f), 2 * FLT_ULP()); 293 test_tol(cexpf, CMPLXL(192, 0x1p-149), 294 CMPLXL(INFINITY, 3.396809344e+38f), 2 * FLT_ULP()); 295 test_tol(cexpf, CMPLXL(120, 0x1.234568p-120), 296 CMPLXL(INFINITY, 1.1163382522e+16f), 2 * FLT_ULP()); 297 test_tol(cexpf, CMPLXL(170, 0x1.234568p-120), 298 CMPLXL(INFINITY, 5.7878851079e+37f), 2 * FLT_ULP()); 299} 300 301ATF_TP_ADD_TCS(tp) 302{ 303 ATF_TP_ADD_TC(tp, zero); 304 ATF_TP_ADD_TC(tp, nan); 305 ATF_TP_ADD_TC(tp, inf); 306 ATF_TP_ADD_TC(tp, reals); 307 ATF_TP_ADD_TC(tp, imaginaries); 308 ATF_TP_ADD_TC(tp, small); 309 ATF_TP_ADD_TC(tp, large); 310 311 return (atf_no_error()); 312} 313