From e44725a8dd99f82f94f37ecff5c0e710c4dba97e Mon Sep 17 00:00:00 2001 From: kettenis Date: Thu, 14 Oct 2021 21:30:00 +0000 Subject: [PATCH] Replace lrint(3), lrintf(3), llrint(3) and llrintf(3) implementations with the same implementation that we're already using for lrintl(3) and llrintl(3). The old implementations were derived from code from NetBSD that didn't pass the lib/libm/msun/lrint_test regress test. NetBSD replaced their implementation with the FreeBSD implementation of this code which we were already using for lrintl(3) and llrintl(3). This fixes the regress test. ok bluhm@, millert@ --- lib/libm/src/s_llrint.c | 14 +++--- lib/libm/src/s_llrintf.c | 14 +++--- lib/libm/src/s_lrint.c | 99 ++++++++++++---------------------------- lib/libm/src/s_lrintf.c | 90 ++++++++++++------------------------ lib/libm/src/s_lrintl.c | 6 +-- 5 files changed, 72 insertions(+), 151 deletions(-) diff --git a/lib/libm/src/s_llrint.c b/lib/libm/src/s_llrint.c index a8d74065b71..85176640490 100644 --- a/lib/libm/src/s_llrint.c +++ b/lib/libm/src/s_llrint.c @@ -1,14 +1,12 @@ -/* $OpenBSD: s_llrint.c,v 1.6 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: llrint.c,v 1.2 2004/10/13 15:18:32 drochner Exp $ */ +/* $OpenBSD: s_llrint.c,v 1.7 2021/10/14 21:30:00 kettenis Exp $ */ /* - * Written by Matthias Drochner . - * Public domain. + * Written by Martynas Venckus. Public domain */ -#define LRINTNAME llrint -#define RESTYPE long long int -#define RESTYPE_MIN LLONG_MIN -#define RESTYPE_MAX LLONG_MAX +#define type double +#define roundit rint +#define dtype long long +#define fn llrint #include "s_lrint.c" diff --git a/lib/libm/src/s_llrintf.c b/lib/libm/src/s_llrintf.c index fdf763693d6..af69ba3fe48 100644 --- a/lib/libm/src/s_llrintf.c +++ b/lib/libm/src/s_llrintf.c @@ -1,14 +1,12 @@ -/* $OpenBSD: s_llrintf.c,v 1.2 2006/09/25 22:16:48 kettenis Exp $ */ -/* $NetBSD: llrintf.c,v 1.2 2004/10/13 15:18:32 drochner Exp $ */ +/* $OpenBSD: s_llrintf.c,v 1.3 2021/10/14 21:30:00 kettenis Exp $ */ /* - * Written by Matthias Drochner . - * Public domain. + * Written by Martynas Venckus. Public domain */ -#define LRINTNAME llrintf -#define RESTYPE long long int -#define RESTYPE_MIN LLONG_MIN -#define RESTYPE_MAX LLONG_MAX +#define type float +#define roundit rintf +#define dtype long long +#define fn llrintf #include "s_lrintf.c" diff --git a/lib/libm/src/s_lrint.c b/lib/libm/src/s_lrint.c index 865cb0e1b40..542ff1da7b7 100644 --- a/lib/libm/src/s_lrint.c +++ b/lib/libm/src/s_lrint.c @@ -1,9 +1,8 @@ -/* $OpenBSD: s_lrint.c,v 1.11 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: lrint.c,v 1.3 2004/10/13 15:18:32 drochner Exp $ */ +/* $OpenBSD: s_lrint.c,v 1.12 2021/10/14 21:30:00 kettenis Exp $ */ /*- - * Copyright (c) 2004 - * Matthias Drochner. All rights reserved. + * Copyright (c) 2005 David Schultz + * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions @@ -27,75 +26,35 @@ * SUCH DAMAGE. */ -#include -#include -#include +#include #include -#include -#include -#include "math_private.h" - -#ifndef LRINTNAME -#define LRINTNAME lrint -#define RESTYPE long int -#define RESTYPE_MIN LONG_MIN -#define RESTYPE_MAX LONG_MAX +#ifndef type +#define type double +#define roundit rint +#define dtype long +#define fn lrint #endif -#define RESTYPE_BITS (sizeof(RESTYPE) * 8) - -static const double -TWO52[2]={ - 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ - -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ -}; - -RESTYPE -LRINTNAME(double x) +/* + * C99 says we should not raise a spurious inexact exception when an + * invalid exception is raised. Unfortunately, the set of inputs + * that overflows depends on the rounding mode when 'dtype' has more + * significant bits than 'type'. Hence, we bend over backwards for the + * sake of correctness; an MD implementation could be more efficient. + */ +dtype +fn(type x) { - u_int32_t i0, i1; - int e, s, shift; - RESTYPE res; - - GET_HIGH_WORD(i0, x); - e = i0 >> DBL_FRACHBITS; - s = e >> DBL_EXPBITS; - e = (e & 0x7ff) - DBL_EXP_BIAS; - - /* 1.0 x 2^31 (or 2^63) is already too large */ - if (e >= (int)RESTYPE_BITS - 1) - return (s ? RESTYPE_MIN : RESTYPE_MAX); /* ??? unspecified */ - - /* >= 2^52 is already an exact integer */ - if (e < DBL_FRACBITS) { - volatile double t = x; /* clip extra precision */ - /* round, using current direction */ - t += TWO52[s]; - t -= TWO52[s]; - x = t; - } - - EXTRACT_WORDS(i0, i1, x); - e = ((i0 >> DBL_FRACHBITS) & 0x7ff) - DBL_EXP_BIAS; - i0 &= 0xfffff; - i0 |= (1 << DBL_FRACHBITS); - - if (e < 0) - return (0); - - shift = e - DBL_FRACBITS; - if (shift >=0) - res = (shift < RESTYPE_BITS ? (RESTYPE)i1 << shift : 0); - else - res = (shift > -RESTYPE_BITS ? (RESTYPE)i1 >> -shift : 0); - shift += 32; - if (shift >=0) - res |= (shift < RESTYPE_BITS ? (RESTYPE)i0 << shift : 0); - else - res |= (shift > -RESTYPE_BITS ? (RESTYPE)i0 >> -shift : 0); - - return (s ? -res : res); + fenv_t env; + dtype d; + + feholdexcept(&env); + d = (dtype)roundit(x); + if (fetestexcept(FE_INVALID)) + feclearexcept(FE_INEXACT); + feupdateenv(&env); + return (d); } -DEF_STD(LRINTNAME); -LDBL_MAYBE_CLONE(LRINTNAME); +DEF_STD(fn); +LDBL_MAYBE_CLONE(fn); diff --git a/lib/libm/src/s_lrintf.c b/lib/libm/src/s_lrintf.c index fae445483e6..a57ee8370bc 100644 --- a/lib/libm/src/s_lrintf.c +++ b/lib/libm/src/s_lrintf.c @@ -1,9 +1,8 @@ -/* $OpenBSD: s_lrintf.c,v 1.6 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: lrintf.c,v 1.3 2004/10/13 15:18:32 drochner Exp $ */ +/* $OpenBSD: s_lrintf.c,v 1.7 2021/10/14 21:30:00 kettenis Exp $ */ /*- - * Copyright (c) 2004 - * Matthias Drochner. All rights reserved. + * Copyright (c) 2005 David Schultz + * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions @@ -27,67 +26,34 @@ * SUCH DAMAGE. */ -#include -#include +#include #include -#include -#include -#include "math_private.h" -#ifndef LRINTNAME -#define LRINTNAME lrintf -#define RESTYPE long int -#define RESTYPE_MIN LONG_MIN -#define RESTYPE_MAX LONG_MAX +#ifndef type +#define type float +#define roundit rintf +#define dtype long +#define fn lrintf #endif -#define RESTYPE_BITS (sizeof(RESTYPE) * 8) - -static const float -TWO23[2]={ - 8.3886080000e+06, /* 0x4b000000 */ - -8.3886080000e+06, /* 0xcb000000 */ -}; - -RESTYPE -LRINTNAME(float x) +/* + * C99 says we should not raise a spurious inexact exception when an + * invalid exception is raised. Unfortunately, the set of inputs + * that overflows depends on the rounding mode when 'dtype' has more + * significant bits than 'type'. Hence, we bend over backwards for the + * sake of correctness; an MD implementation could be more efficient. + */ +dtype +fn(type x) { - u_int32_t i0; - int e, s, shift; - RESTYPE res; - - GET_FLOAT_WORD(i0, x); - e = i0 >> SNG_FRACBITS; - s = e >> SNG_EXPBITS; - e = (e & 0xff) - SNG_EXP_BIAS; - - /* 1.0 x 2^31 (or 2^63) is already too large */ - if (e >= (int)RESTYPE_BITS - 1) - return (s ? RESTYPE_MIN : RESTYPE_MAX); /* ??? unspecified */ - - /* >= 2^23 is already an exact integer */ - if (e < SNG_FRACBITS) { - volatile float t = x; /* clip extra precision */ - /* round, using current direction */ - t += TWO23[s]; - t -= TWO23[s]; - x = t; - } - - GET_FLOAT_WORD(i0, x); - e = ((i0 >> SNG_FRACBITS) & 0xff) - SNG_EXP_BIAS; - i0 &= 0x7fffff; - i0 |= (1 << SNG_FRACBITS); - - if (e < 0) - return (0); - - shift = e - SNG_FRACBITS; - if (shift >=0) - res = (shift < RESTYPE_BITS ? (RESTYPE)i0 << shift : 0); - else - res = (shift > -RESTYPE_BITS ? (RESTYPE)i0 >> -shift : 0); - - return (s ? -res : res); + fenv_t env; + dtype d; + + feholdexcept(&env); + d = (dtype)roundit(x); + if (fetestexcept(FE_INVALID)) + feclearexcept(FE_INEXACT); + feupdateenv(&env); + return (d); } -DEF_STD(LRINTNAME); +DEF_STD(fn); diff --git a/lib/libm/src/s_lrintl.c b/lib/libm/src/s_lrintl.c index 62b3155c2b5..d317b30fc14 100644 --- a/lib/libm/src/s_lrintl.c +++ b/lib/libm/src/s_lrintl.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_lrintl.c,v 1.3 2019/03/15 05:42:38 kevlo Exp $ */ +/* $OpenBSD: s_lrintl.c,v 1.4 2021/10/14 21:30:00 kettenis Exp $ */ /*- * Copyright (c) 2005 David Schultz @@ -31,9 +31,9 @@ #ifndef type #define type long double -#define roundit rintl +#define roundit rintl #define dtype long -#define fn lrintl +#define fn lrintl #endif /* -- 2.20.1