244 lines
6.3 KiB
C
244 lines
6.3 KiB
C
#ifdef HAVE_CONFIG_H
|
|
#include <config.h>
|
|
#endif
|
|
|
|
#include <math.h>
|
|
|
|
/*
|
|
* Integral value nearest to x in the current rounding direction (C23
|
|
* 7.12.9.4), all three precisions. rint differs from nearbyint only in
|
|
* that it may raise the inexact exception; never raising it is also
|
|
* conforming, and the const attribute promises no errno path.
|
|
*
|
|
* The __builtin_rint forms are tempting (they fold to in-line code on
|
|
* this target), but GCC warns -Winfinite-recursion when the enclosing
|
|
* function carries the same name as the library symbol the builtin would
|
|
* fall back to (a function named rint whose body is __builtin_rint),
|
|
* which the -Wall -Wextra -pedantic build gate forbids. Each function is
|
|
* therefore implemented directly on the IEEE 754 bit pattern with
|
|
* round-to-nearest-even, exactly like nearbyint (see nearbyint.c): no
|
|
* <fenv.h> exists in vlibc yet, so the hardware default round-to-nearest-
|
|
* even is the only reachable rounding mode. rint(±0) is ±0, rint(-0.5)
|
|
* is -0.0 (tie to even zero), and ±Inf/NaN pass through unchanged.
|
|
*/
|
|
|
|
/*
|
|
* The double format: bit 63 is the sign, bits 62..52 the exponent biased
|
|
* by 1023, bits 51..0 the fraction. For e in [1023, 1074] the low
|
|
* (1075 - e) bits of the fraction word are the fractional part; the
|
|
* significand is 53 bits wide, so in the e == 1023 binade the tie-even
|
|
* test looks at the implicit bit (the only half-way value there, 1.5,
|
|
* rounds up to 2).
|
|
*/
|
|
static double
|
|
rint_d(double x)
|
|
{
|
|
unsigned long long bits;
|
|
unsigned long long kept;
|
|
unsigned long long frac;
|
|
unsigned long long half;
|
|
int e;
|
|
int shift;
|
|
|
|
__builtin_memcpy(&bits, &x, sizeof bits);
|
|
e = (int)((bits >> 52) & 0x7ff);
|
|
|
|
if (e >= 1075)
|
|
{
|
|
return x;
|
|
}
|
|
if (e < 1023)
|
|
{
|
|
if ((bits & ~(1ULL << 63)) == 0)
|
|
{
|
|
return x;
|
|
}
|
|
if (e == 1022)
|
|
{
|
|
/* [1/2, 1): 0.5 itself is a tie toward even zero; anything
|
|
* above it rounds to ±1. */
|
|
if ((bits & 0xFFFFFFFFFFFFFULL) == 0)
|
|
{
|
|
bits &= 1ULL << 63;
|
|
__builtin_memcpy(&x, &bits, sizeof x);
|
|
return x;
|
|
}
|
|
bits = (bits & (1ULL << 63)) | 0x3FF0000000000000ULL;
|
|
__builtin_memcpy(&x, &bits, sizeof x);
|
|
return x;
|
|
}
|
|
bits &= 1ULL << 63;
|
|
__builtin_memcpy(&x, &bits, sizeof x);
|
|
return x;
|
|
}
|
|
|
|
shift = 1075 - e;
|
|
half = 1ULL << (shift - 1);
|
|
frac = bits & ((1ULL << shift) - 1ULL);
|
|
kept = bits & ~((1ULL << shift) - 1ULL);
|
|
if (frac > half || (frac == half && (shift == 52 || ((kept >> shift) & 1ULL) != 0)))
|
|
{
|
|
kept += 1ULL << shift;
|
|
}
|
|
__builtin_memcpy(&x, &kept, sizeof x);
|
|
return x;
|
|
}
|
|
|
|
/*
|
|
* The float format: bit 31 is the sign, bits 30..23 the exponent biased
|
|
* by 127, bits 22..0 the fraction; the significand is 24 bits wide. For
|
|
* e in [127, 149] the low (150 - e) fraction bits are fractional, and in
|
|
* the e == 127 binade the tie-even test looks at the implicit bit.
|
|
*/
|
|
static float
|
|
rint_f(float x)
|
|
{
|
|
unsigned int bits;
|
|
unsigned int kept;
|
|
unsigned int frac;
|
|
unsigned int half;
|
|
int e;
|
|
int shift;
|
|
|
|
__builtin_memcpy(&bits, &x, sizeof bits);
|
|
e = (int)((bits >> 23) & 0xff);
|
|
|
|
if (e >= 150)
|
|
{
|
|
return x;
|
|
}
|
|
if (e < 127)
|
|
{
|
|
if ((bits & ~(1U << 31)) == 0)
|
|
{
|
|
return x;
|
|
}
|
|
if (e == 126)
|
|
{
|
|
if ((bits & 0x7FFFFFU) == 0)
|
|
{
|
|
bits &= 1U << 31;
|
|
__builtin_memcpy(&x, &bits, sizeof x);
|
|
return x;
|
|
}
|
|
bits = (bits & (1U << 31)) | 0x3F800000U;
|
|
__builtin_memcpy(&x, &bits, sizeof x);
|
|
return x;
|
|
}
|
|
bits &= 1U << 31;
|
|
__builtin_memcpy(&x, &bits, sizeof x);
|
|
return x;
|
|
}
|
|
|
|
shift = 150 - e;
|
|
half = 1U << (shift - 1);
|
|
frac = bits & ((1U << shift) - 1U);
|
|
kept = bits & ~((1U << shift) - 1U);
|
|
if (frac > half || (frac == half && (shift == 23 || ((kept >> shift) & 1U) != 0)))
|
|
{
|
|
kept += 1U << shift;
|
|
}
|
|
__builtin_memcpy(&x, &kept, sizeof x);
|
|
return x;
|
|
}
|
|
|
|
/*
|
|
* The x86 80-bit extended format: 64 significand bits m (the integer bit
|
|
* is explicit, so the tie-even test is always m's kept LSB) and a
|
|
* sign/exponent word se. For e in [16383, 16445] the low (16446 - e)
|
|
* bits of m are fractional. The step may overflow m when the kept
|
|
* significand is all ones; the carry then moves the value to the next
|
|
* binade.
|
|
*/
|
|
static long double
|
|
rint_ld(long double x)
|
|
{
|
|
struct
|
|
{
|
|
unsigned long long m;
|
|
unsigned short se;
|
|
} p;
|
|
unsigned long long frac;
|
|
unsigned long long half;
|
|
int e;
|
|
int shift;
|
|
|
|
__builtin_memcpy(&p, &x, sizeof p);
|
|
e = p.se & 0x7fff;
|
|
|
|
if (e >= 16446)
|
|
{
|
|
return x;
|
|
}
|
|
if (e < 16383)
|
|
{
|
|
if (p.m == 0)
|
|
{
|
|
return x;
|
|
}
|
|
if (e == 16382)
|
|
{
|
|
if (p.m == 0x8000000000000000ULL)
|
|
{
|
|
/* Exactly 0.5: tie toward even zero. */
|
|
p.m = 0;
|
|
p.se &= 0x8000;
|
|
__builtin_memcpy(&x, &p, sizeof p);
|
|
return x;
|
|
}
|
|
p.m = 0x8000000000000000ULL;
|
|
p.se = (p.se & 0x8000) | 16383;
|
|
__builtin_memcpy(&x, &p, sizeof p);
|
|
return x;
|
|
}
|
|
p.m = 0;
|
|
p.se &= 0x8000;
|
|
__builtin_memcpy(&x, &p, sizeof p);
|
|
return x;
|
|
}
|
|
|
|
shift = 16446 - e;
|
|
half = 1ULL << (shift - 1);
|
|
frac = p.m & ((1ULL << shift) - 1ULL);
|
|
p.m &= ~((1ULL << shift) - 1ULL);
|
|
if (frac > half || (frac == half && ((p.m >> shift) & 1ULL) != 0))
|
|
{
|
|
p.m += 1ULL << shift;
|
|
if (p.m == 0)
|
|
{
|
|
/* Kept significand was all ones: carry to the next binade. */
|
|
p.m = 0x8000000000000000ULL;
|
|
p.se = (p.se & 0x8000) | (unsigned short)(e + 1);
|
|
}
|
|
}
|
|
__builtin_memcpy(&x, &p, sizeof p);
|
|
return x;
|
|
}
|
|
|
|
/*
|
|
* As rint, for a float argument.
|
|
*/
|
|
float
|
|
rintf(float x)
|
|
{
|
|
return rint_f(x);
|
|
}
|
|
|
|
/*
|
|
* As rint, for a double argument.
|
|
*/
|
|
double
|
|
rint(double x)
|
|
{
|
|
return rint_d(x);
|
|
}
|
|
|
|
/*
|
|
* As rint, for a long double argument.
|
|
*/
|
|
long double
|
|
rintl(long double x)
|
|
{
|
|
return rint_ld(x);
|
|
}
|