mirror of
https://https.git.savannah.gnu.org/git/gnulib.git
synced 2026-06-15 15:25:49 +00:00
Collin Funk
118 lines
3.6 KiB
C
118 lines
3.6 KiB
C
/* Round according to the current rounding mode.
|
|
Copyright (C) 2007, 2010-2026 Free Software Foundation, Inc.
|
|
|
|
This file is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU Lesser General Public License as
|
|
published by the Free Software Foundation, either version 3 of the
|
|
License, or (at your option) any later version.
|
|
|
|
This file is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU Lesser General Public License for more details.
|
|
|
|
You should have received a copy of the GNU Lesser General Public License
|
|
along with this program. If not, see <https://www.gnu.org/licenses/>. */
|
|
|
|
#if ! defined USE_LONG_DOUBLE
|
|
# include <config.h>
|
|
#endif
|
|
|
|
/* Specification. */
|
|
#include <math.h>
|
|
|
|
#include <float.h>
|
|
#include <stdlib.h>
|
|
|
|
#undef MIN
|
|
|
|
#ifdef USE_LONG_DOUBLE
|
|
# define RINT rintl
|
|
# define DOUBLE long double
|
|
# define MANT_DIG LDBL_MANT_DIG
|
|
# define MIN LDBL_MIN
|
|
# define L_(literal) literal##L
|
|
#elif ! defined USE_FLOAT
|
|
# define RINT rint
|
|
# define DOUBLE double
|
|
# define MANT_DIG DBL_MANT_DIG
|
|
# define MIN DBL_MIN
|
|
# define L_(literal) literal
|
|
#else /* defined USE_FLOAT */
|
|
# define RINT rintf
|
|
# define DOUBLE float
|
|
# define MANT_DIG FLT_MANT_DIG
|
|
# define MIN FLT_MIN
|
|
# define L_(literal) literal##f
|
|
#endif
|
|
|
|
/* -0.0. See minus-zero.h. */
|
|
#if defined __hpux || defined __ICC
|
|
# define MINUS_ZERO (-MIN * MIN)
|
|
#else
|
|
# define MINUS_ZERO L_(-0.0)
|
|
#endif
|
|
|
|
/* MSVC with option -fp:strict refuses to compile constant initializers that
|
|
contain floating-point operations. Pacify this compiler. */
|
|
#if defined _MSC_VER && !defined __clang__
|
|
# pragma fenv_access (off)
|
|
#endif
|
|
|
|
DOUBLE
|
|
RINT (DOUBLE x)
|
|
{
|
|
/* 2^(MANT_DIG-1). */
|
|
static const DOUBLE TWO_MANT_DIG =
|
|
/* Assume MANT_DIG <= 5 * 31.
|
|
Use the identity
|
|
n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5). */
|
|
(DOUBLE) (1U << ((MANT_DIG - 1) / 5))
|
|
* (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5))
|
|
* (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5))
|
|
* (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5))
|
|
* (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5));
|
|
|
|
/* The use of 'volatile' guarantees that excess precision bits are dropped at
|
|
each addition step and before the following comparison at the caller's
|
|
site. It is necessary on x86 systems where double-floats are not IEEE
|
|
compliant by default, to avoid that the results become platform and
|
|
compiler option dependent. 'volatile' is a portable alternative to gcc's
|
|
-ffloat-store option. */
|
|
volatile DOUBLE z = x;
|
|
|
|
/* Consider the current rounding mode, cf.
|
|
<https://pubs.opengroup.org/onlinepubs/9699919799/basedefs/float.h.html>.
|
|
This implementation supports only rounds-to-nearest. */
|
|
if (FLT_ROUNDS != 1)
|
|
abort ();
|
|
|
|
if (z > L_(0.0))
|
|
{
|
|
/* Avoid rounding error for x = 0.5 - 2^(-MANT_DIG-1). */
|
|
if (z < L_(0.5))
|
|
z = L_(0.0);
|
|
/* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */
|
|
else if (z < TWO_MANT_DIG)
|
|
{
|
|
/* Round to the next integer. */
|
|
z += TWO_MANT_DIG;
|
|
z -= TWO_MANT_DIG;
|
|
}
|
|
}
|
|
else if (z < L_(0.0))
|
|
{
|
|
/* Avoid rounding error for x = -(0.5 - 2^(-MANT_DIG-1)). */
|
|
if (z > - L_(0.5))
|
|
z = MINUS_ZERO;
|
|
/* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */
|
|
else if (z > -TWO_MANT_DIG)
|
|
{
|
|
/* Round to the next integer. */
|
|
z -= TWO_MANT_DIG;
|
|
z += TWO_MANT_DIG;
|
|
}
|
|
}
|
|
return z;
|
|
}
|