t_hypot.c revision 1.4
1/* $NetBSD: t_hypot.c,v 1.4 2024/05/11 20:09:13 riastradh Exp $ */ 2 3/*- 4 * Copyright (c) 2016 The NetBSD Foundation, Inc. 5 * All rights reserved. 6 * 7 * Redistribution and use in source and binary forms, with or without 8 * modification, are permitted provided that the following conditions 9 * are met: 10 * 1. Redistributions of source code must retain the above copyright 11 * notice, this list of conditions and the following disclaimer. 12 * 2. Redistributions in binary form must reproduce the above copyright 13 * notice, this list of conditions and the following disclaimer in the 14 * documentation and/or other materials provided with the distribution. 15 * 16 * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS 17 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 18 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 19 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS 20 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 21 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 22 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 23 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 24 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 25 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 26 * POSSIBILITY OF SUCH DAMAGE. 27 */ 28 29#include <atf-c.h> 30#include <float.h> 31#include <math.h> 32 33#define CHECK_EQ(i, hypot, a, b, c) \ 34 ATF_CHECK_MSG(hypot(a, b) == (c), \ 35 "[%u] %s(%a, %a)=%a, expected %a", \ 36 (i), #hypot, (a), (b), hypot(a, b), (c)) 37 38#define CHECKL_EQ(i, hypot, a, b, c) \ 39 ATF_CHECK_MSG(hypot(a, b) == (c), \ 40 "[%u] %s(%La, %La)=%La, expected %La", \ 41 (i), #hypot, (a), (b), hypot(a, b), (c)) 42 43static const float trivial_casesf[] = { 44 0, 45#ifdef __FLT_HAS_DENORM__ 46 __FLT_DENORM_MIN__, 47 2*__FLT_DENORM_MIN__, 48 3*__FLT_DENORM_MIN__, 49 FLT_MIN - 3*__FLT_DENORM_MIN__, 50 FLT_MIN - 2*__FLT_DENORM_MIN__, 51 FLT_MIN - __FLT_DENORM_MIN__, 52#endif 53 FLT_MIN, 54 FLT_MIN*(1 + FLT_EPSILON), 55 FLT_MIN*(1 + 2*FLT_EPSILON), 56 2*FLT_MIN, 57 FLT_EPSILON/2, 58 FLT_EPSILON, 59 2*FLT_EPSILON, 60 1 - 3*FLT_EPSILON/2, 61 1 - 2*FLT_EPSILON/2, 62 1 - FLT_EPSILON/2, 63 1, 64 1 + FLT_EPSILON, 65 1 + 2*FLT_EPSILON, 66 1 + 3*FLT_EPSILON, 67 1.5 - 3*FLT_EPSILON, 68 1.5 - 2*FLT_EPSILON, 69 1.5 - FLT_EPSILON, 70 1.5, 71 1.5 + FLT_EPSILON, 72 1.5 + 2*FLT_EPSILON, 73 1.5 + 3*FLT_EPSILON, 74 2, 75 0.5/FLT_EPSILON - 0.5, 76 0.5/FLT_EPSILON, 77 0.5/FLT_EPSILON + 0.5, 78 1/FLT_EPSILON, 79 FLT_MAX, 80 INFINITY, 81}; 82 83static const double trivial_cases[] = { 84 0, 85#ifdef __DBL_HAS_DENORM__ 86 __DBL_DENORM_MIN__, 87 2*__DBL_DENORM_MIN__, 88 3*__DBL_DENORM_MIN__, 89 DBL_MIN - 3*__DBL_DENORM_MIN__, 90 DBL_MIN - 2*__DBL_DENORM_MIN__, 91 DBL_MIN - __DBL_DENORM_MIN__, 92#endif 93 DBL_MIN, 94 DBL_MIN*(1 + DBL_EPSILON), 95 DBL_MIN*(1 + 2*DBL_EPSILON), 96 2*DBL_MIN, 97 DBL_EPSILON/2, 98 DBL_EPSILON, 99 2*DBL_EPSILON, 100 1 - 3*DBL_EPSILON/2, 101 1 - 2*DBL_EPSILON/2, 102 1 - DBL_EPSILON/2, 103 1, 104 1 + DBL_EPSILON, 105 1 + 2*DBL_EPSILON, 106 1 + 3*DBL_EPSILON, 107 1.5 - 3*DBL_EPSILON, 108 1.5 - 2*DBL_EPSILON, 109 1.5 - DBL_EPSILON, 110 1.5, 111 1.5 + DBL_EPSILON, 112 1.5 + 2*DBL_EPSILON, 113 1.5 + 3*DBL_EPSILON, 114 2, 115 1/FLT_EPSILON - 0.5, 116 1/FLT_EPSILON, 117 1/FLT_EPSILON + 0.5, 118 0.5/DBL_EPSILON - 0.5, 119 0.5/DBL_EPSILON, 120 0.5/DBL_EPSILON + 0.5, 121 1/DBL_EPSILON, 122 DBL_MAX, 123 INFINITY, 124}; 125 126static const long double trivial_casesl[] = { 127 0, 128#ifdef __LDBL_HAS_DENORM__ 129 __LDBL_DENORM_MIN__, 130 2*__LDBL_DENORM_MIN__, 131 3*__LDBL_DENORM_MIN__, 132 LDBL_MIN - 3*__LDBL_DENORM_MIN__, 133 LDBL_MIN - 2*__LDBL_DENORM_MIN__, 134 LDBL_MIN - __LDBL_DENORM_MIN__, 135#endif 136 LDBL_MIN, 137 LDBL_MIN*(1 + LDBL_EPSILON), 138 LDBL_MIN*(1 + 2*LDBL_EPSILON), 139 2*LDBL_MIN, 140 LDBL_EPSILON/2, 141 LDBL_EPSILON, 142 2*LDBL_EPSILON, 143 1 - 3*LDBL_EPSILON/2, 144 1 - 2*LDBL_EPSILON/2, 145 1 - LDBL_EPSILON/2, 146 1, 147 1 + LDBL_EPSILON, 148 1 + 2*LDBL_EPSILON, 149 1 + 3*LDBL_EPSILON, 150 1.5 - 3*LDBL_EPSILON, 151 1.5 - 2*LDBL_EPSILON, 152 1.5 - LDBL_EPSILON, 153 1.5, 154 1.5 + LDBL_EPSILON, 155 1.5 + 2*LDBL_EPSILON, 156 1.5 + 3*LDBL_EPSILON, 157 2, 158 1/FLT_EPSILON - 0.5, 159 1/FLT_EPSILON, 160 1/FLT_EPSILON + 0.5, 161#ifdef __HAVE_LONG_DOUBLE 162 1/DBL_EPSILON - 0.5L, 163 1/DBL_EPSILON, 164 1/DBL_EPSILON + 0.5L, 165#endif 166 0.5/LDBL_EPSILON - 0.5, 167 0.5/LDBL_EPSILON, 168 0.5/LDBL_EPSILON + 0.5, 169 1/LDBL_EPSILON, 170 LDBL_MAX, 171 INFINITY, 172}; 173 174ATF_TC(hypotf_trivial); 175ATF_TC_HEAD(hypotf_trivial, tc) 176{ 177 atf_tc_set_md_var(tc, "descr", "hypotf(x,0) and hypotf(0,x)"); 178} 179ATF_TC_BODY(hypotf_trivial, tc) 180{ 181 unsigned i; 182 183 for (i = 0; i < __arraycount(trivial_casesf); i++) { 184 volatile float x = trivial_casesf[i]; 185 186 CHECK_EQ(i, hypotf, x, +0., x); 187 CHECK_EQ(i, hypotf, x, -0., x); 188 CHECK_EQ(i, hypotf, +0., x, x); 189 CHECK_EQ(i, hypotf, -0., x, x); 190 CHECK_EQ(i, hypotf, -x, +0., x); 191 CHECK_EQ(i, hypotf, -x, -0., x); 192 CHECK_EQ(i, hypotf, +0., -x, x); 193 CHECK_EQ(i, hypotf, -0., -x, x); 194 } 195} 196 197ATF_TC(hypot_trivial); 198ATF_TC_HEAD(hypot_trivial, tc) 199{ 200 atf_tc_set_md_var(tc, "descr", "hypot(x,0) and hypot(0,x)"); 201} 202ATF_TC_BODY(hypot_trivial, tc) 203{ 204 unsigned i; 205 206 for (i = 0; i < __arraycount(trivial_casesf); i++) { 207 volatile double x = trivial_casesf[i]; 208 209 CHECK_EQ(i, hypot, x, +0., x); 210 CHECK_EQ(i, hypot, x, -0., x); 211 CHECK_EQ(i, hypot, +0., x, x); 212 CHECK_EQ(i, hypot, -0., x, x); 213 CHECK_EQ(i, hypot, -x, +0., x); 214 CHECK_EQ(i, hypot, -x, -0., x); 215 CHECK_EQ(i, hypot, +0., -x, x); 216 CHECK_EQ(i, hypot, -0., -x, x); 217 } 218 219 for (i = 0; i < __arraycount(trivial_cases); i++) { 220 volatile double x = trivial_cases[i]; 221 222 CHECK_EQ(i, hypot, x, +0., x); 223 CHECK_EQ(i, hypot, x, -0., x); 224 CHECK_EQ(i, hypot, +0., x, x); 225 CHECK_EQ(i, hypot, -0., x, x); 226 CHECK_EQ(i, hypot, -x, +0., x); 227 CHECK_EQ(i, hypot, -x, -0., x); 228 CHECK_EQ(i, hypot, +0., -x, x); 229 CHECK_EQ(i, hypot, -0., -x, x); 230 } 231} 232 233ATF_TC(hypotl_trivial); 234ATF_TC_HEAD(hypotl_trivial, tc) 235{ 236 atf_tc_set_md_var(tc, "descr", "hypotl(x,0) and hypotl(0,x)"); 237} 238ATF_TC_BODY(hypotl_trivial, tc) 239{ 240 unsigned i; 241 242 for (i = 0; i < __arraycount(trivial_casesf); i++) { 243 volatile long double x = trivial_casesf[i]; 244 245 CHECKL_EQ(i, hypotl, x, +0.L, x); 246 CHECKL_EQ(i, hypotl, x, -0.L, x); 247 CHECKL_EQ(i, hypotl, +0.L, x, x); 248 CHECKL_EQ(i, hypotl, -0.L, x, x); 249 CHECKL_EQ(i, hypotl, -x, +0.L, x); 250 CHECKL_EQ(i, hypotl, -x, -0.L, x); 251 CHECKL_EQ(i, hypotl, +0.L, -x, x); 252 CHECKL_EQ(i, hypotl, -0.L, -x, x); 253 } 254 255 for (i = 0; i < __arraycount(trivial_cases); i++) { 256 volatile long double x = trivial_cases[i]; 257 258 CHECKL_EQ(i, hypotl, x, +0.L, x); 259 CHECKL_EQ(i, hypotl, x, -0.L, x); 260 CHECKL_EQ(i, hypotl, +0.L, x, x); 261 CHECKL_EQ(i, hypotl, -0.L, x, x); 262 CHECKL_EQ(i, hypotl, -x, +0.L, x); 263 CHECKL_EQ(i, hypotl, -x, -0.L, x); 264 CHECKL_EQ(i, hypotl, +0.L, -x, x); 265 CHECKL_EQ(i, hypotl, -0.L, -x, x); 266 } 267 268#if __HAVE_LONG_DOUBLE + 0 == 128 269 atf_tc_expect_fail("PR lib/58245: hypotl is broken on ld128 ports"); 270#endif 271 272 for (i = 0; i < __arraycount(trivial_casesl); i++) { 273 volatile long double x = trivial_casesl[i]; 274 275 CHECKL_EQ(i, hypotl, x, +0.L, x); 276 CHECKL_EQ(i, hypotl, x, -0.L, x); 277 CHECKL_EQ(i, hypotl, +0.L, x, x); 278 CHECKL_EQ(i, hypotl, -0.L, x, x); 279 CHECKL_EQ(i, hypotl, -x, +0.L, x); 280 CHECKL_EQ(i, hypotl, -x, -0.L, x); 281 CHECKL_EQ(i, hypotl, +0.L, -x, x); 282 CHECKL_EQ(i, hypotl, -0.L, -x, x); 283 } 284} 285 286__CTASSERT(FLT_MANT_DIG >= 24); 287static const struct { 288 float a, b, c; 289} exact_casesf[] = { 290 { 3, 4, 5 }, 291 { 5, 12, 13 }, 292 { 9, 12, 15 }, 293 { 0x1001, 0x801000, 0x801001 }, 294 { 4248257, 1130976, 4396225 }, 295}; 296 297__CTASSERT(DBL_MANT_DIG >= 53); 298static const struct { 299 double a, b, c; 300} exact_cases[] = { 301 { 3378249543467007, 4505248894795776, 5631148868747265 }, 302 { 0x7ffffff, 0x1ffffff8000000, 0x1ffffff8000001 }, 303#if DBL_MANT_DIG >= 56 304 { 13514123525517439, 18018830919909120, 22523538851237761 }, 305 { 0x1fffffff, 0x1ffffffe0000000, 0x1ffffffe0000001 }, 306#endif 307}; 308 309#if LDBL_MANT_DIG >= 64 310static const struct { 311 long double a, b, c; 312} exact_casesl[] = { 313 { 3458976450080784639, 4611968592949214720, 5764960744407842561 }, 314 { 0x1ffffffff, 0x1fffffffe00000000p0L, 0x1fffffffe00000001p0L }, 315#if LDBL_MANT_DIG >= 113 316 { 973555668229277869436257492279295.L, 317 1298074224305703705819019479072768.L, 318 1622592780382129686316970078625793.L }, 319 { 0x1ffffffffffffff, 320 0x1fffffffffffffe00000000000000p0L, 321 0x1fffffffffffffe00000000000001p0L }, 322#endif 323}; 324#endif 325 326ATF_TC(hypotf_exact); 327ATF_TC_HEAD(hypotf_exact, tc) 328{ 329 atf_tc_set_md_var(tc, "descr", "hypotf on scaled Pythagorean triples"); 330} 331ATF_TC_BODY(hypotf_exact, tc) 332{ 333 unsigned i; 334 335 for (i = 0; i < __arraycount(exact_casesf); i++) { 336 int s; 337 338 for (s = FLT_MIN_EXP; 339 s < FLT_MAX_EXP - FLT_MANT_DIG; 340 s += (FLT_MAX_EXP - FLT_MANT_DIG - FLT_MIN_EXP)/5) { 341 volatile double a = ldexpf(exact_casesf[i].a, s); 342 volatile double b = ldexpf(exact_casesf[i].b, s); 343 float c = ldexpf(exact_casesf[i].c, s); 344 345 CHECK_EQ(i, hypot, a, b, c); 346 CHECK_EQ(i, hypot, b, a, c); 347 CHECK_EQ(i, hypot, a, -b, c); 348 CHECK_EQ(i, hypot, b, -a, c); 349 CHECK_EQ(i, hypot, -a, b, c); 350 CHECK_EQ(i, hypot, -b, a, c); 351 CHECK_EQ(i, hypot, -a, -b, c); 352 CHECK_EQ(i, hypot, -b, -a, c); 353 } 354 } 355} 356 357ATF_TC(hypot_exact); 358ATF_TC_HEAD(hypot_exact, tc) 359{ 360 atf_tc_set_md_var(tc, "descr", "hypot on scaled Pythagorean triples"); 361} 362ATF_TC_BODY(hypot_exact, tc) 363{ 364 unsigned i; 365 366 for (i = 0; i < __arraycount(exact_casesf); i++) { 367 int s; 368 369 for (s = DBL_MIN_EXP; 370 s < DBL_MAX_EXP - DBL_MANT_DIG; 371 s += (DBL_MAX_EXP - DBL_MANT_DIG - DBL_MIN_EXP)/5) { 372 volatile double a = ldexp(exact_casesf[i].a, s); 373 volatile double b = ldexp(exact_casesf[i].b, s); 374 double c = ldexp(exact_casesf[i].c, s); 375 376 CHECK_EQ(i, hypot, a, b, c); 377 CHECK_EQ(i, hypot, b, a, c); 378 CHECK_EQ(i, hypot, a, -b, c); 379 CHECK_EQ(i, hypot, b, -a, c); 380 CHECK_EQ(i, hypot, -a, b, c); 381 CHECK_EQ(i, hypot, -b, a, c); 382 CHECK_EQ(i, hypot, -a, -b, c); 383 CHECK_EQ(i, hypot, -b, -a, c); 384 } 385 } 386 387 for (i = 0; i < __arraycount(exact_cases); i++) { 388 int s; 389 390 for (s = DBL_MIN_EXP; 391 s < DBL_MAX_EXP - DBL_MANT_DIG; 392 s += (DBL_MAX_EXP - DBL_MANT_DIG - DBL_MIN_EXP)/5) { 393 volatile double a = ldexp(exact_cases[i].a, s); 394 volatile double b = ldexp(exact_cases[i].b, s); 395 double c = ldexp(exact_cases[i].c, s); 396 397 CHECK_EQ(i, hypot, a, b, c); 398 CHECK_EQ(i, hypot, b, a, c); 399 CHECK_EQ(i, hypot, a, -b, c); 400 CHECK_EQ(i, hypot, b, -a, c); 401 CHECK_EQ(i, hypot, -a, b, c); 402 CHECK_EQ(i, hypot, -b, a, c); 403 CHECK_EQ(i, hypot, -a, -b, c); 404 CHECK_EQ(i, hypot, -b, -a, c); 405 } 406 } 407} 408 409ATF_TC(hypotl_exact); 410ATF_TC_HEAD(hypotl_exact, tc) 411{ 412 atf_tc_set_md_var(tc, "descr", "hypotl on scaled Pythagorean triples"); 413} 414ATF_TC_BODY(hypotl_exact, tc) 415{ 416 unsigned i; 417 418 for (i = 0; i < __arraycount(exact_casesf); i++) { 419 int s; 420 421 for (s = LDBL_MIN_EXP; 422 s < LDBL_MAX_EXP - LDBL_MANT_DIG; 423 s += (LDBL_MAX_EXP - LDBL_MANT_DIG - LDBL_MIN_EXP)/5) { 424 volatile long double a = ldexpl(exact_casesf[i].a, s); 425 volatile long double b = ldexpl(exact_casesf[i].b, s); 426 long double c = ldexpl(exact_casesf[i].c, s); 427 428 CHECKL_EQ(i, hypotl, a, b, c); 429 CHECKL_EQ(i, hypotl, b, a, c); 430 CHECKL_EQ(i, hypotl, a, -b, c); 431 CHECKL_EQ(i, hypotl, b, -a, c); 432 CHECKL_EQ(i, hypotl, -a, b, c); 433 CHECKL_EQ(i, hypotl, -b, a, c); 434 CHECKL_EQ(i, hypotl, -a, -b, c); 435 CHECKL_EQ(i, hypotl, -b, -a, c); 436 } 437 } 438 439 for (i = 0; i < __arraycount(exact_cases); i++) { 440 int s; 441 442 for (s = LDBL_MIN_EXP; 443 s < LDBL_MAX_EXP - LDBL_MANT_DIG; 444 s += (LDBL_MAX_EXP - LDBL_MANT_DIG - LDBL_MIN_EXP)/5) { 445 volatile long double a = ldexpl(exact_cases[i].a, s); 446 volatile long double b = ldexpl(exact_cases[i].b, s); 447 long double c = ldexpl(exact_cases[i].c, s); 448 449 CHECKL_EQ(i, hypotl, a, b, c); 450 CHECKL_EQ(i, hypotl, b, a, c); 451 CHECKL_EQ(i, hypotl, a, -b, c); 452 CHECKL_EQ(i, hypotl, b, -a, c); 453 CHECKL_EQ(i, hypotl, -a, b, c); 454 CHECKL_EQ(i, hypotl, -b, a, c); 455 CHECKL_EQ(i, hypotl, -a, -b, c); 456 CHECKL_EQ(i, hypotl, -b, -a, c); 457 } 458 } 459 460#if __HAVE_LONG_DOUBLE + 0 == 128 461 atf_tc_expect_fail("PR lib/58245: hypotl is broken on ld128 ports"); 462#endif 463 464#if LDBL_MANT_DIG >= 64 465 for (i = 0; i < __arraycount(exact_casesl); i++) { 466 int s; 467 468 for (s = LDBL_MIN_EXP; 469 s < LDBL_MAX_EXP - LDBL_MANT_DIG; 470 s += (LDBL_MAX_EXP - LDBL_MANT_DIG - LDBL_MIN_EXP)/5) { 471 volatile long double a = ldexpl(exact_casesl[i].a, s); 472 volatile long double b = ldexpl(exact_casesl[i].b, s); 473 long double c = ldexpl(exact_casesl[i].c, s); 474 475 CHECKL_EQ(i, hypotl, a, b, c); 476 CHECKL_EQ(i, hypotl, b, a, c); 477 CHECKL_EQ(i, hypotl, a, -b, c); 478 CHECKL_EQ(i, hypotl, b, -a, c); 479 CHECKL_EQ(i, hypotl, -a, b, c); 480 CHECKL_EQ(i, hypotl, -b, a, c); 481 CHECKL_EQ(i, hypotl, -a, -b, c); 482 CHECKL_EQ(i, hypotl, -b, -a, c); 483 } 484 } 485#endif 486} 487 488ATF_TC(pr50698); 489ATF_TC_HEAD(pr50698, tc) 490{ 491 atf_tc_set_md_var(tc, "descr", "Check for the bug of PR lib/50698"); 492} 493 494ATF_TC_BODY(pr50698, tc) 495{ 496 volatile float a = 1e-18f; 497 float val = hypotf(a, a); 498 ATF_CHECK(!isinf(val)); 499 ATF_CHECK(!isnan(val)); 500} 501 502ATF_TP_ADD_TCS(tp) 503{ 504 505 ATF_TP_ADD_TC(tp, hypot_exact); 506 ATF_TP_ADD_TC(tp, hypot_trivial); 507 ATF_TP_ADD_TC(tp, hypotf_exact); 508 ATF_TP_ADD_TC(tp, hypotf_trivial); 509 ATF_TP_ADD_TC(tp, hypotl_exact); 510 ATF_TP_ADD_TC(tp, hypotl_trivial); 511 ATF_TP_ADD_TC(tp, pr50698); 512 513 return atf_no_error(); 514} 515