From b731a02db45f80b10f35227459c3190bd34ef7d7 Mon Sep 17 00:00:00 2001 From: huntedbytheirs Date: Sun, 6 Sep 2026 01:07:31 -0400 Subject: [PATCH] feat(math): abs/round/trunc/frexp/ldexp/scalbn/copysign/fmin/fmax/fmod --- include/math.h | 470 ++++++++++++++++++++++++++++++++++++++++ include/tgmath.h | 71 +++++- src/math/ceil.c | 173 +++++++++++++++ src/math/copysign.c | 74 +++++++ src/math/fabs.c | 41 ++++ src/math/fdim.c | 71 ++++++ src/math/floor.c | 207 ++++++++++++++++++ src/math/fmax.c | 261 ++++++++++++++++++++++ src/math/fmin.c | 267 +++++++++++++++++++++++ src/math/fmod.c | 423 ++++++++++++++++++++++++++++++++++++ src/math/frexp.c | 173 +++++++++++++++ src/math/ilogb.c | 126 +++++++++++ src/math/ldexp.c | 76 +++++++ src/math/llrint.c | 45 ++++ src/math/llround.c | 73 +++++++ src/math/logb.c | 123 +++++++++++ src/math/lrint.c | 49 +++++ src/math/lround.c | 64 ++++++ src/math/math_impl.h | 341 +++++++++++++++++++++++++++++ src/math/modf.c | 215 ++++++++++++++++++ src/math/nearbyint.c | 248 +++++++++++++++++++++ src/math/remainder.c | 284 ++++++++++++++++++++++++ src/math/remquo.c | 364 +++++++++++++++++++++++++++++++ src/math/rint.c | 243 +++++++++++++++++++++ src/math/round.c | 226 +++++++++++++++++++ src/math/scalbn.c | 121 +++++++++++ src/math/trunc.c | 127 +++++++++++ tests/test_math_basic.c | 342 +++++++++++++++++++++++++++++ tests/test_math_s1.c | 264 ++++++++++++++++++++++ tests/test_math_s2.c | 357 ++++++++++++++++++++++++++++++ 30 files changed, 5912 insertions(+), 7 deletions(-) create mode 100644 include/math.h create mode 100644 src/math/ceil.c create mode 100644 src/math/copysign.c create mode 100644 src/math/fabs.c create mode 100644 src/math/fdim.c create mode 100644 src/math/floor.c create mode 100644 src/math/fmax.c create mode 100644 src/math/fmin.c create mode 100644 src/math/fmod.c create mode 100644 src/math/frexp.c create mode 100644 src/math/ilogb.c create mode 100644 src/math/ldexp.c create mode 100644 src/math/llrint.c create mode 100644 src/math/llround.c create mode 100644 src/math/logb.c create mode 100644 src/math/lrint.c create mode 100644 src/math/lround.c create mode 100644 src/math/math_impl.h create mode 100644 src/math/modf.c create mode 100644 src/math/nearbyint.c create mode 100644 src/math/remainder.c create mode 100644 src/math/remquo.c create mode 100644 src/math/rint.c create mode 100644 src/math/round.c create mode 100644 src/math/scalbn.c create mode 100644 src/math/trunc.c create mode 100644 tests/test_math_basic.c create mode 100644 tests/test_math_s1.c create mode 100644 tests/test_math_s2.c diff --git a/include/math.h b/include/math.h new file mode 100644 index 0000000..395a776 --- /dev/null +++ b/include/math.h @@ -0,0 +1,470 @@ +#ifndef VLIBC_MATH_H +#define VLIBC_MATH_H + +/* + * vlibc — . + * + * This header is ISO C core and is present in every profile. It carries the + * whole todo-39 arithmetic inventory: the classification macros + * (fpclassify/isnan/isinf/isfinite/isnormal/signbit), the constants + * (HUGE_VAL/INFINITY/NAN, FP_*, math_errhandling), and the 25 basic + * real-function families (fabs, copysign, floor, ceil, trunc, round, rint, + * nearbyint, lrint, llrint, frexp, ldexp, modf, scalbn, scalbln, fmin, + * fmax, fdim, fmod, remainder, remquo, ilogb, logb, lround, llround), each + * spelled in all three precisions (float `...f`, double, long double + * `...l`). + * + * Every declaration below lands NOW so that the header is stable across the + * todo-39 implementation slices; the slices (wip commits) fill in the + * src/math/ definitions family by family and never touch this header again. + * Later math todos extend this header in place the same way stdlib.h is + * extended: todo 40 adds the exp/log/pow families, todo 41 the trig and + * hyperbolic families, todo 42 erf/erfc/lgamma/tgamma/fma/nextafter/nan, + * and todo 43 provides separately. + * + * Edge-case behavior (IEEE 754-2008): floor/ceil/trunc round exact results + * with correct signed-zero and NaN/Inf handling; fabs/copysign/fmin/fmax/ + * fdim never raise an exception on their own and preserve NaN payloads + * where the standard allows; the classification macros are + * exception-free. Functions that per POSIX can set errno (ldexp, scalbn, + * scalbln, fmod, remainder, remquo, and the fraction-splitters frexp/modf + * with their pointer outputs) carry no const attribute below so the + * compiler never hoists or elides an errno-setting call. Functions that + * are pure (never set errno in their defined domain — fabs, copysign, + * floor, ceil, trunc, round, rint, nearbyint, fmin, fmax, fdim, ilogb, + * logb) are declared const and fold away in static links. + */ + +#include + +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +/* ---- Constants ---- */ + +/* + * Positive infinity, as double/float/long double (C23 7.12.1p4). HUGE_VAL + * is also the "overflowed" return value of the strto* conversions. + */ +#define HUGE_VAL __builtin_huge_val() +#define HUGE_VALF __builtin_huge_valf() +#define HUGE_VALL __builtin_huge_vall() + +/* Positive infinity as a float (C23 7.12.1p5). */ +#define INFINITY __builtin_inff() + +/* A quiet NaN as a float (C23 7.12.1p6); the payload is implementation + * defined but always a NaN of the quiet kind. */ +#define NAN __builtin_nanf("") + +/* + * Return values of fpclassify (C23 7.12.3.1): the numbers need only be + * distinct positive values; the classification macros below hand these + * exact constants to __builtin_fpclassify so the two always agree. + */ +#define FP_NAN 0 +#define FP_INFINITE 1 +#define FP_ZERO 2 +#define FP_SUBNORMAL 3 +#define FP_NORMAL 4 + +/* + * Return values of ilogb (C23 7.12.6.5p4): FP_ILOGB0 for a zero argument + * and FP_ILOGBNAN for a NaN or infinite argument, each the indicated + * sentinel. + */ +#define FP_ILOGB0 INT_MIN +#define FP_ILOGBNAN INT_MAX + +/* + * Which error mechanisms the library reports through (C23 7.12.1p3): + * MATH_ERRNO and MATH_ERREXCEPT are both in effect (the domain/range + * errno values are set where POSIX requires and the corresponding + * floating-point exceptions are raised by the hardware). + */ +#define MATH_ERRNO 1 +#define MATH_ERREXCEPT 2 +#define math_errhandling 3 + +/* ---- Classification macros (type-generic, exception-free) ---- */ + +/* + * Classify x as NaN, infinite, zero, subnormal, or normal, returning the + * matching FP_* constant. __builtin_fpclassify is a GCC type-generic + * builtin that evaluates its floating argument exactly once and never + * traps, so no _Generic dispatch is needed here. + */ +#define fpclassify(x) \ + __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, (x)) + +/* + * The predicate macros below are thin wrappers over the corresponding + * GCC type-generic builtins. Each builtin evaluates its argument once and + * returns an int; integer arguments classify as finite, non-zero, and + * normal (they convert exactly), so isnan(3) is 0 and isinf(3) is 0. + */ +#define isnan(x) __builtin_isnan(x) +#define isinf(x) __builtin_isinf(x) +#define isfinite(x) __builtin_isfinite(x) +#define isnormal(x) __builtin_isnormal(x) +#define signbit(x) __builtin_signbit(x) + +/* ---- Basic arithmetic families (todo 39) ---- */ + +/* + * Absolute value of x (C23 7.12.7.2). fabs(±0) is +0, fabs(±Inf) is +Inf, + * and fabs(NaN) is a NaN. Pure: no domain, no exception. + */ +__attribute__((const)) float +fabsf(float x); + +__attribute__((const)) double +fabs(double x); + +__attribute__((const)) long double +fabsl(long double x); + +/* + * A value with the magnitude of x and the sign of y (C23 7.12.7.3). + * copysign(±0, y) carries y's sign; a NaN x keeps its payload but takes + * y's sign bit. Pure. + */ +__attribute__((const)) float +copysignf(float x, float y); + +__attribute__((const)) double +copysign(double x, double y); + +__attribute__((const)) long double +copysignl(long double x, long double y); + +/* + * The largest integral value not greater than x (C23 7.12.9.2). + * floor(±0) is ±0, floor(-0.5) is -1.0, floor(±Inf) and floor(NaN) + * return their argument unchanged. Pure and exact (no inexact + * exception), so it folds. + */ +__attribute__((const)) float +floorf(float x); + +__attribute__((const)) double +floor(double x); + +__attribute__((const)) long double +floorl(long double x); + +/* + * The smallest integral value not less than x (C23 7.12.9.1). + * ceil(±0) is ±0, ceil(-0.5) is -0.0 (a real negative zero), ceil(±Inf) + * and ceil(NaN) return their argument unchanged. Pure and exact. + */ +__attribute__((const)) float +ceilf(float x); + +__attribute__((const)) double +ceil(double x); + +__attribute__((const)) long double +ceill(long double x); + +/* + * The integral value nearest to x in the direction of zero (C23 7.12.9.3). + * trunc(-0.7) is -0.0, trunc(0.7) is 0.0, and ±Inf/NaN pass through. + * Pure and exact. + */ +__attribute__((const)) float +truncf(float x); + +__attribute__((const)) double +trunc(double x); + +__attribute__((const)) long double +truncl(long double x); + +/* + * The integral value nearest to x, with halfway cases rounded away from + * zero (C23 7.12.9.6). round(±0) is ±0, round(-0.5) is -1.0, and ±Inf/ + * NaN pass through. Pure: round never raises the inexact exception. + */ +__attribute__((const)) float +roundf(float x); + +__attribute__((const)) double +round(double x); + +__attribute__((const)) long double +roundl(long double x); + +/* + * The integral value nearest to x in the current rounding direction + * (C23 7.12.9.4). rint may raise the inexact exception; it is still pure + * in the const sense because it never sets errno and reads no memory. + */ +__attribute__((const)) float +rintf(float x); + +__attribute__((const)) double +rint(double x); + +__attribute__((const)) long double +rintl(long double x); + +/* + * As rint, but guaranteed never to raise the inexact exception (C23 + * 7.12.9.5). Pure and exact. + */ +__attribute__((const)) float +nearbyintf(float x); + +__attribute__((const)) double +nearbyint(double x); + +__attribute__((const)) long double +nearbyintl(long double x); + +/* + * The nearest integral value to x in the current rounding direction, + * returned as long (C23 7.12.9.7). A result outside the range of long is + * a range error (the return value is unspecified and errno may be set), + * so no const attribute. + */ +long +lrintf(float x); + +long +lrint(double x); + +long +lrintl(long double x); + +/* + * As lrint, returned as long long (C23 7.12.9.8). Range errors as for + * lrint; no const attribute. + */ +long long +llrintf(float x); + +long long +llrint(double x); + +long long +llrintl(long double x); + +/* + * Split x into a fraction f in [1/2, 1) (or 0) and an integer exponent + * *exp such that x == f * 2^(*exp) (C23 7.12.6.4). Zero returns ±0 with + * *exp 0; Inf/NaN return x with an unspecified *exp. Writes *exp, so no + * const attribute. + */ +float +frexpf(float x, int *exp); + +double +frexp(double x, int *exp); + +long double +frexpl(long double x, int *exp); + +/* + * x times 2^n (C23 7.12.6.3): the inverse of frexp. A result too large + * to represent is a range error returning ±HUGE_VAL with errno ERANGE; + * no const attribute. + */ +float +ldexpf(float x, int n); + +double +ldexp(double x, int n); + +long double +ldexpl(long double x, int n); + +/* + * Split x into an integral part stored in *iptr and a fractional part + * returned (C23 7.12.6.5); both have x's sign, so modf(-1.5, &i) puts + * -1.0 in i and returns -0.5. Writes *iptr, so no const attribute. + */ +float +modff(float x, float *iptr); + +double +modf(double x, double *iptr); + +long double +modfl(long double x, long double *iptr); + +/* + * x * FLT_RADIX^n with FLT_RADIX 2 (C23 7.12.6.6): the scalb* functions + * differ from ldexp only in the exponent argument's type. A result too + * large to represent is a range error returning ±HUGE_VAL* with errno + * ERANGE; no const attribute. + */ +float +scalbnf(float x, int n); + +double +scalbn(double x, int n); + +long double +scalbnl(long double x, int n); + +float +scalblnf(float x, long n); + +double +scalbln(double x, long n); + +long double +scalblnl(long double x, long n); + +/* + * The smaller of x and y (C23 7.12.12.4), returning -0.0 when the + * arguments are +0.0 and -0.0. A NaN argument is ignored in favor of the + * numeric one; two NaNs return a NaN. Pure. + */ +__attribute__((const)) float +fminf(float x, float y); + +__attribute__((const)) double +fmin(double x, double y); + +__attribute__((const)) long double +fminl(long double x, long double y); + +/* + * The larger of x and y (C23 7.12.12.3), returning +0.0 when the + * arguments are +0.0 and -0.0. NaN handling as fmin. Pure. + */ +__attribute__((const)) float +fmaxf(float x, float y); + +__attribute__((const)) double +fmax(double x, double y); + +__attribute__((const)) long double +fmaxl(long double x, long double y); + +/* + * The positive difference x - y when x > y and +0.0 otherwise (C23 + * 7.12.12.2); fdim(x, NaN) and fdim(NaN, x) return a NaN. Pure in the + * domain-error-free sense (an overflow may still raise an exception + * through the hardware, which const does not model). + */ +__attribute__((const)) float +fdimf(float x, float y); + +__attribute__((const)) double +fdim(double x, double y); + +__attribute__((const)) long double +fdiml(long double x, long double y); + +/* + * The floating-point remainder x - n*y, where n is x/y truncated toward + * zero (C23 7.12.10.1); the result therefore has x's sign. fmod(x, ±0) is + * a domain error returning NaN with errno EDOM; ±Inf/x is likewise a + * domain error. errno can be set, so no const attribute. + */ +float +fmodf(float x, float y); + +double +fmod(double x, double y); + +long double +fmodl(long double x, long double y); + +/* + * The IEEE remainder x - n*y, where n is x/y rounded to the nearest + * integer (ties to even) (C23 7.12.10.2); |result| <= |y|/2. Domain + * errors as fmod; errno can be set, so no const attribute. + */ +float +remainderf(float x, float y); + +double +remainder(double x, double y); + +long double +remainderl(long double x, long double y); + +/* + * As remainder, additionally storing the low bits of the integer + * quotient n in *quo (C23 7.12.10.3). Writes *quo, so no const + * attribute. + */ +float +remquof(float x, float y, int *quo); + +double +remquo(double x, double y, int *quo); + +long double +remquol(long double x, long double y, int *quo); + +/* + * The signed exponent of x as an int (C23 7.12.6.5p1): ilogb(x) is + * floor(log2 |x|) for a nonzero finite x. ilogb(±0) returns FP_ILOGB0, + * ilogb(±Inf) and ilogb(NaN) return FP_ILOGBNAN. Pure: the sentinels are + * returned without touching errno. + */ +__attribute__((const)) int +ilogbf(float x); + +__attribute__((const)) int +ilogb(double x); + +__attribute__((const)) int +ilogbl(long double x); + +/* + * The signed exponent of x as a floating-point value (C23 7.12.6.6): + * logb(±0) is -Inf and logb(±Inf) is +Inf, each raising the division- + * by-zero/invalid exception through the hardware but without an errno + * path in the representable domain. Pure. + */ +__attribute__((const)) float +logbf(float x); + +__attribute__((const)) double +logb(double x); + +__attribute__((const)) long double +logbl(long double x); + +/* + * The nearest integral value to x, with halfway cases rounded away from + * zero, returned as long (C23 7.12.9.9). A result outside the range of + * long is a range error; no const attribute. + */ +long +lroundf(float x); + +long +lround(double x); + +long +lroundl(long double x); + +/* + * As lround, returned as long long (C23 7.12.9.10). Range errors as for + * lround; no const attribute. + */ +long long +llroundf(float x); + +long long +llround(double x); + +long long +llroundl(long double x); + +#ifdef __cplusplus +} +#endif + +#endif /* VLIBC_MATH_H */ diff --git a/include/tgmath.h b/include/tgmath.h index 09744f6..6b09e49 100644 --- a/include/tgmath.h +++ b/include/tgmath.h @@ -4,15 +4,72 @@ /* * vlibc — . * - * Type-generic math. STUB: the type-generic dispatch over is - * IMPLEMENTED by the first math todo (#39), not here. This stub exists so - * that including alongside the other headers compiles cleanly - * today; it deliberately defines no type-generic macros until #39 lands - * the function inventory it dispatches to. + * Type-generic math (C23 7.25). Each macro inspects the type of its + * controlling expression and dispatches to the matching precision of the + * real function: a float argument selects the `...f` variant, a + * long double argument the `...l` variant, and everything else (double, + * any integer type, and the _Complex types, whose imaginary part is + * dropped by the usual argument conversion) the unsuffixed double + * variant. The controlling expression is the first floating argument, so + * the two-argument forms dispatch on x and require y to be convertible to + * x's type. * - * This header is ISO C core and is present in every profile. + * This header is ISO C core and is present in every profile. The + * classification macros (fpclassify/isnan/isinf/isfinite/isnormal/signbit) + * need no dispatch of their own — already defines them over + * GCC's type-generic builtins, so including below makes them + * work unchanged for every real and complex type. */ -#include +#include + +/* fabs */ +#define fabs(x) _Generic((x), float: fabsf, long double: fabsl, default: fabs)(x) + +/* copysign */ +#define copysign(x, y) \ + _Generic((x), float: copysignf, long double: copysignl, default: copysign)(x, y) + +/* rounding toward zero / ±Inf / current mode (floor, ceil, trunc, round, + * rint, nearbyint) */ +#define floor(x) _Generic((x), float: floorf, long double: floorl, default: floor)(x) +#define ceil(x) _Generic((x), float: ceilf, long double: ceill, default: ceil)(x) +#define trunc(x) _Generic((x), float: truncf, long double: truncl, default: trunc)(x) +#define round(x) _Generic((x), float: roundf, long double: roundl, default: round)(x) +#define rint(x) _Generic((x), float: rintf, long double: rintl, default: rint)(x) +#define nearbyint(x) \ + _Generic((x), float: nearbyintf, long double: nearbyintl, default: nearbyint)(x) + +/* integer-returning rounding (lrint, llrint, lround, llround) */ +#define lrint(x) _Generic((x), float: lrintf, long double: lrintl, default: lrint)(x) +#define llrint(x) _Generic((x), float: llrintf, long double: llrintl, default: llrint)(x) +#define lround(x) _Generic((x), float: lroundf, long double: lroundl, default: lround)(x) +#define llround(x) _Generic((x), float: llroundf, long double: llroundl, default: llround)(x) + +/* fraction/exponent splitting (frexp, modf) and exponent scaling (ldexp, + * scalbn, scalbln): the second argument is an int/long or a pointer to + * the first argument's type and follows x's precision automatically */ +#define frexp(x, e) _Generic((x), float: frexpf, long double: frexpl, default: frexp)(x, e) +#define modf(x, i) _Generic((x), float: modff, long double: modfl, default: modf)(x, i) +#define ldexp(x, n) _Generic((x), float: ldexpf, long double: ldexpl, default: ldexp)(x, n) +#define scalbn(x, n) _Generic((x), float: scalbnf, long double: scalbnl, default: scalbn)(x, n) +#define scalbln(x, n) _Generic((x), float: scalblnf, long double: scalblnl, default: scalbln)(x, n) + +/* minimum/maximum/difference */ +#define fmin(x, y) _Generic((x), float: fminf, long double: fminl, default: fmin)(x, y) +#define fmax(x, y) _Generic((x), float: fmaxf, long double: fmaxl, default: fmax)(x, y) +#define fdim(x, y) _Generic((x), float: fdimf, long double: fdiml, default: fdim)(x, y) + +/* remainder (fmod, remainder, remquo: the trailing int* follows x's + * precision) */ +#define fmod(x, y) _Generic((x), float: fmodf, long double: fmodl, default: fmod)(x, y) +#define remainder(x, y) \ + _Generic((x), float: remainderf, long double: remainderl, default: remainder)(x, y) +#define remquo(x, y, q) \ + _Generic((x), float: remquof, long double: remquol, default: remquo)(x, y, q) + +/* exponent reading (ilogb returns int, logb returns x's precision) */ +#define ilogb(x) _Generic((x), float: ilogbf, long double: ilogbl, default: ilogb)(x) +#define logb(x) _Generic((x), float: logbf, long double: logbl, default: logb)(x) #endif /* VLIBC_TGMATH_H */ diff --git a/src/math/ceil.c b/src/math/ceil.c new file mode 100644 index 0000000..79f4478 --- /dev/null +++ b/src/math/ceil.c @@ -0,0 +1,173 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Smallest integral value not less than x (C23 7.12.9.1), all three + * precisions. + * + * As with floor (see floor.c), GCC does not fold the __builtin_ceil forms + * on this target, so each function works on the IEEE 754 bit pattern. + * Clearing the fraction bits rounds toward zero, which is already the + * correct direction for a negative argument; a positive argument with a + * nonzero fraction must then step up by one (an exact integer + 1.0). + * ceil(-0.3) is a real -0.0, and ±0/±Inf/NaN pass through unchanged. + */ + +static double +ceil_d(double x) +{ + const unsigned long long sign_mask = 1ULL << 63; + unsigned long long bits; + unsigned long long frac; + 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 & ~sign_mask) == 0) + { + return x; + } + return (bits & sign_mask) != 0 ? -0.0 : 1.0; + } + + shift = 1075 - e; + frac = (1ULL << shift) - 1; + if ((bits & sign_mask) != 0) + { + bits &= ~frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + if ((bits & frac) == 0) + { + return x; + } + bits &= ~frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x + 1.0; +} + +static float +ceil_f(float x) +{ + const unsigned int sign_mask = 1U << 31; + unsigned int bits; + unsigned int frac; + 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 & ~sign_mask) == 0) + { + return x; + } + return (bits & sign_mask) != 0 ? -0.0f : 1.0f; + } + + shift = 150 - e; + frac = (1U << shift) - 1U; + if ((bits & sign_mask) != 0) + { + bits &= ~frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + if ((bits & frac) == 0) + { + return x; + } + bits &= ~frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x + 1.0f; +} + +static long double +ceil_ld(long double x) +{ + struct + { + unsigned long long m; + unsigned short se; + } p; + unsigned long long frac; + 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 && (p.se & 0x7fff) == 0) + { + return x; + } + return (p.se & 0x8000) != 0 ? -0.0L : 1.0L; + } + + shift = 63 - (e - 16383); + frac = (1ULL << shift) - 1ULL; + if ((p.se & 0x8000) != 0) + { + p.m &= ~frac; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + if ((p.m & frac) == 0) + { + return x; + } + p.m &= ~frac; + __builtin_memcpy(&x, &p, sizeof p); + return x + 1.0L; +} + +/* + * As ceil, for a float argument. + */ +float +ceilf(float x) +{ + return ceil_f(x); +} + +/* + * As ceil, for a double argument. + */ +double +ceil(double x) +{ + return ceil_d(x); +} + +/* + * As ceil, for a long double argument. + */ +long double +ceill(long double x) +{ + return ceil_ld(x); +} diff --git a/src/math/copysign.c b/src/math/copysign.c new file mode 100644 index 0000000..6bafdd8 --- /dev/null +++ b/src/math/copysign.c @@ -0,0 +1,74 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Copy the sign of y onto the magnitude of x, all three precisions (C23 + * 7.12.7.3). Each function below splices the sign bit of y into the IEEE + * 754 bit pattern of x: the magnitude of x is untouched (so a NaN keeps + * its payload) and copysign(±0, y) carries y's sign. GCC's + * __builtin_copysign forms do fold to andp/orp pairs on this target, but + * the compiler diagnoses the __builtin_ call inside the identically named + * function as infinite recursion, so the splice is written out directly; + * it compiles to the same two instructions. + */ + +/* + * A float with the magnitude of x and the sign of y. + */ +float +copysignf(float x, float y) // NOLINT(bugprone-easily-swappable-parameters) +{ + unsigned int xb; + unsigned int yb; + + __builtin_memcpy(&xb, &x, sizeof xb); + __builtin_memcpy(&yb, &y, sizeof yb); + xb = (xb & 0x7fffffffU) | (yb & 0x80000000U); + __builtin_memcpy(&x, &xb, sizeof xb); + return x; +} + +/* + * A double with the magnitude of x and the sign of y. + */ +double +copysign(double x, double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + unsigned long long xb; + unsigned long long yb; + + __builtin_memcpy(&xb, &x, sizeof xb); + __builtin_memcpy(&yb, &y, sizeof yb); + xb = (xb & 0x7fffffffffffffffULL) | (yb & 0x8000000000000000ULL); + __builtin_memcpy(&x, &xb, sizeof xb); + return x; +} + +/* + * A long double with the magnitude of x and the sign of y. The x86 80-bit + * extended format keeps the sign in bit 15 of the sign/exponent word at + * bytes 8..9, so only that bit is spliced. + */ +long double +copysignl(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + struct + { + unsigned long long m; + unsigned short se; + } xp; + struct + { + unsigned long long m; + unsigned short se; + } yp; + + __builtin_memcpy(&xp, &x, sizeof xp); + __builtin_memcpy(&yp, &y, sizeof yp); + xp.se = (unsigned short)((xp.se & 0x7fff) | (yp.se & 0x8000)); + __builtin_memcpy(&x, &xp, sizeof xp); + return x; +} diff --git a/src/math/fabs.c b/src/math/fabs.c new file mode 100644 index 0000000..5a8d9c1 --- /dev/null +++ b/src/math/fabs.c @@ -0,0 +1,41 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Absolute value, all three precisions (C23 7.12.7.2). GCC folds + * __builtin_fabs/__builtin_fabsf/__builtin_fabsl into a single + * sign-clearing SSE/x87 instruction at every optimization level (verified + * at -O0 and -O2), so the call never recurses and no errno path exists. + * The sign-bit clear is exact: fabs(-0.0) is +0.0, fabs(±Inf) is +Inf, + * and a NaN keeps its payload. + */ + +/* + * Absolute value of x as a float. + */ +float +fabsf(float x) +{ + return __builtin_fabsf(x); +} + +/* + * Absolute value of x as a double. + */ +double +fabs(double x) +{ + return __builtin_fabs(x); +} + +/* + * Absolute value of x as a long double. + */ +long double +fabsl(long double x) +{ + return __builtin_fabsl(x); +} diff --git a/src/math/fdim.c b/src/math/fdim.c new file mode 100644 index 0000000..11e6dde --- /dev/null +++ b/src/math/fdim.c @@ -0,0 +1,71 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * The positive difference x - y when x > y and +0.0 otherwise (C23 + * 7.12.12.2), all three precisions. + * + * C23 keeps fdim(x, NaN) and fdim(NaN, x) NaN, so the two arguments are + * screened before the ordering test; a NaN result is the quiet NaN that + * the x + y addition below produces (matching the host glibc). With both + * operands numeric, x > y is an ordinary comparison (no exceptions) and + * the single subtraction x - y is the whole computation: an exact + * representable difference stays exact, and a difference too large for the + * format overflows through the hardware into +Inf exactly as glibc's does, + * with the overflow flag raised and no extra help needed. x <= y -- the + * signed-zero and equal cases included -- returns a plain +0.0. + */ + +/* + * fdim of two floats. + */ +float +fdimf(float x, float y) // NOLINT(bugprone-easily-swappable-parameters) +{ + if (isnan(x) || isnan(y)) + { + return x + y; + } + if (x > y) + { + return x - y; + } + return 0.0f; +} + +/* + * fdim of two doubles. + */ +double +fdim(double x, double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + if (isnan(x) || isnan(y)) + { + return x + y; + } + if (x > y) + { + return x - y; + } + return 0.0; +} + +/* + * fdim of two long doubles. + */ +long double +fdiml(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + if (isnan(x) || isnan(y)) + { + return x + y; + } + if (x > y) + { + return x - y; + } + return 0.0L; +} diff --git a/src/math/floor.c b/src/math/floor.c new file mode 100644 index 0000000..01c95fe --- /dev/null +++ b/src/math/floor.c @@ -0,0 +1,207 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Largest integral value not greater than x (C23 7.12.9.2), all three + * precisions. + * + * GCC emits an external floor/floorf/floorl call for the __builtin_ forms + * at every optimization level on this target (no SSE4.1 in the default + * -march, so the roundsd expansion is unavailable and the builtin is not + * folded), so each function below is implemented directly on the IEEE 754 + * bit pattern instead: clearing the fraction bits rounds toward zero, and + * a negative argument with a nonzero fraction must then step down by one. + * The clearing and the exact integer - 1.0 are both exact, so no rounding + * mode and no floating-point exception is involved; ±0/±Inf/NaN pass + * through and signed zero is preserved. + */ + +/* + * The double format: bit 63 is the sign, bits 62..52 the exponent biased + * by 1023, bits 51..0 the fraction. A finite value with exponent e has + * |x| in [2^(e-1023), 2^(e-1022)) and a fraction only when e - 1023 < 52, + * i.e. e < 1075; for e in [1023, 1074] exactly the low (1075 - e) bits of + * the fraction word are the fractional part. + */ +static double +floor_d(double x) +{ + const unsigned long long sign_mask = 1ULL << 63; + unsigned long long bits; + unsigned long long frac; + int e; + int shift; + + __builtin_memcpy(&bits, &x, sizeof bits); + e = (int)((bits >> 52) & 0x7ff); + + /* |x| >= 2^52 is already integral; Inf (e == 0x7ff) and NaN must pass + * through unchanged as well. */ + if (e >= 1075) + { + return x; + } + + /* |x| < 1: floor is +0 for a nonnegative x and -1 for a negative one, + * except that ±0 is its own floor. */ + if (e < 1023) + { + if ((bits & ~sign_mask) == 0) + { + return x; + } + return (bits & sign_mask) != 0 ? -1.0 : 0.0; + } + + shift = 1075 - e; + frac = (1ULL << shift) - 1; + if ((bits & sign_mask) == 0) + { + bits &= ~frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + + /* Negative: round toward zero, then take one more step down when a + * fraction was dropped. The truncation leaves an integer-valued + * double, so subtracting 1.0 is exact. */ + if ((bits & frac) == 0) + { + return x; + } + bits &= ~frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x - 1.0; +} + +/* + * The float format: bit 31 is the sign, bits 30..23 the exponent biased + * by 127, bits 22..0 the fraction. Fraction bits exist exactly when the + * exponent e is in [127, 149]; the low (150 - e) bits are fractional. + */ +static float +floor_f(float x) +{ + const unsigned int sign_mask = 1U << 31; + unsigned int bits; + unsigned int frac; + 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 & ~sign_mask) == 0) + { + return x; + } + return (bits & sign_mask) != 0 ? -1.0f : 0.0f; + } + + shift = 150 - e; + frac = (1U << shift) - 1U; + if ((bits & sign_mask) == 0) + { + bits &= ~frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + if ((bits & frac) == 0) + { + return x; + } + bits &= ~frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x - 1.0f; +} + +/* + * The x86 80-bit extended format: 64 significand bits m (the integer bit + * is explicit) in bytes 0..7 and a sign/exponent word se in bytes 8..9, + * with the sign in bit 15 and the exponent (biased by 16383) in bits + * 14..0. A value m * 2^(e - 16446) has fractional bits only when + * e - 16383 < 63, i.e. e < 16446; for e in [16383, 16445] exactly the low + * (63 - (e - 16383)) bits of m are fractional. + */ +static long double +floor_ld(long double x) +{ + struct + { + unsigned long long m; + unsigned short se; + } p; + unsigned long long frac; + int e; + int shift; + + __builtin_memcpy(&p, &x, sizeof p); + e = p.se & 0x7fff; + + /* |x| >= 2^63 is already integral; Inf (e == 0x7fff) and NaN pass + * through unchanged as well. */ + if (e >= 16446) + { + return x; + } + if (e < 16383) + { + if (p.m == 0 && (p.se & 0x7fff) == 0) + { + return x; + } + return (p.se & 0x8000) != 0 ? -1.0L : 0.0L; + } + + shift = 63 - (e - 16383); + frac = (1ULL << shift) - 1ULL; + if ((p.se & 0x8000) == 0) + { + p.m &= ~frac; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + if ((p.m & frac) == 0) + { + return x; + } + p.m &= ~frac; + __builtin_memcpy(&x, &p, sizeof p); + return x - 1.0L; +} + +/* + * As floor, for a float argument. + */ +float +floorf(float x) +{ + return floor_f(x); +} + +/* + * As floor, for a double argument. + */ +double +floor(double x) +{ + return floor_d(x); +} + +/* + * As floor, for a long double argument. + */ +long double +floorl(long double x) +{ + return floor_ld(x); +} diff --git a/src/math/fmax.c b/src/math/fmax.c new file mode 100644 index 0000000..9f45af2 --- /dev/null +++ b/src/math/fmax.c @@ -0,0 +1,261 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * The larger of x and y (C23 7.12.12.3), all three precisions. + * + * The mirror image of fmin (see fmin.c for the full reasoning): a quiet + * NaN argument is ignored in favor of the numeric one, a signaling NaN + * argument makes the result that signaling NaN quieted, and a +0.0/-0.0 + * pair returns +0.0. Only the equal-argument tie rule (and, with it, the + * direction of the zero preference) differs from fmin; the numeric + * comparisons, the NaN screens and the word-level quiet splice are + * otherwise identical, so the implementation shares the shape of fmin.c + * and is likewise pure. + */ + +/* + * The double format: bit 63 is the sign, bits 62..52 the exponent biased + * by 1023, bits 51..0 the fraction; exponent 0x7ff with a nonzero fraction + * is a NaN, quiet bit 0x0008000000000000. + */ +static int +fmax_nan_d(unsigned long long b) +{ + return ((b >> 52) & 0x7ff) == 0x7ff && (b & 0xFFFFFFFFFFFFFULL) != 0; +} + +static int +fmax_snan_d(unsigned long long b) +{ + return fmax_nan_d(b) && (b & 0x0008000000000000ULL) == 0; +} + +static double +fmax_mm_d(double x, double y) +{ + unsigned long long xb; + unsigned long long yb; + int xneg; + int yneg; + + __builtin_memcpy(&xb, &x, sizeof xb); + __builtin_memcpy(&yb, &y, sizeof yb); + + if (fmax_snan_d(xb)) + { + xb |= 0x0008000000000000ULL; + __builtin_memcpy(&x, &xb, sizeof xb); + return x; + } + if (fmax_snan_d(yb)) + { + yb |= 0x0008000000000000ULL; + __builtin_memcpy(&y, &yb, sizeof yb); + return y; + } + if (fmax_nan_d(xb)) + { + return fmax_nan_d(yb) ? x : y; + } + if (fmax_nan_d(yb)) + { + return x; + } + + if (x < y) + { + return y; + } + if (y < x) + { + return x; + } + + /* Equal, so possibly a +0/-0 pair: fmax prefers +0. */ + if (x == 0.0) + { + xneg = (int)(xb >> 63); + yneg = (int)(yb >> 63); + if (xneg != yneg) + { + return xneg ? y : x; + } + } + return x; +} + +/* + * The float format: bit 31 the sign, bits 30..23 the exponent biased by + * 127, bits 22..0 the fraction; quiet bit 0x00400000. + */ +static int +fmax_nan_f(unsigned int b) +{ + return ((b >> 23) & 0xff) == 0xff && (b & 0x7FFFFFU) != 0; +} + +static int +fmax_snan_f(unsigned int b) +{ + return fmax_nan_f(b) && (b & 0x00400000U) == 0; +} + +static float +fmax_mm_f(float x, float y) +{ + unsigned int xb; + unsigned int yb; + int xneg; + int yneg; + + __builtin_memcpy(&xb, &x, sizeof xb); + __builtin_memcpy(&yb, &y, sizeof yb); + + if (fmax_snan_f(xb)) + { + xb |= 0x00400000U; + __builtin_memcpy(&x, &xb, sizeof xb); + return x; + } + if (fmax_snan_f(yb)) + { + yb |= 0x00400000U; + __builtin_memcpy(&y, &yb, sizeof yb); + return y; + } + if (fmax_nan_f(xb)) + { + return fmax_nan_f(yb) ? x : y; + } + if (fmax_nan_f(yb)) + { + return x; + } + + if (x < y) + { + return y; + } + if (y < x) + { + return x; + } + if (x == 0.0f) + { + xneg = (int)(xb >> 31); + yneg = (int)(yb >> 31); + if (xneg != yneg) + { + return xneg ? y : x; + } + } + return x; +} + +/* + * The x86 80-bit extended format: 64-bit significand m (explicit integer + * bit) in bytes 0..7, sign/exponent word se in bytes 8..9 with the sign in + * bit 15 and the exponent biased by 16383 in bits 14..0. NaN has + * (se & 0x7fff) == 0x7fff with a nonzero fraction; the quiet bit is bit 62 + * of m. + */ +struct fmax_ld_word +{ + unsigned long long m; + unsigned short se; +}; + +static int +fmax_nan_l(struct fmax_ld_word p) +{ + return (p.se & 0x7fff) == 0x7fff && (p.m & 0x7FFFFFFFFFFFFFFFULL) != 0; +} + +static int +fmax_snan_l(struct fmax_ld_word p) +{ + return fmax_nan_l(p) && (p.m & (1ULL << 62)) == 0; +} + +static long double +fmax_mm_l(long double x, long double y) +{ + struct fmax_ld_word p; + struct fmax_ld_word q; + int xneg; + int yneg; + + __builtin_memcpy(&p, &x, sizeof p); + __builtin_memcpy(&q, &y, sizeof q); + + if (fmax_snan_l(p)) + { + p.m |= 1ULL << 62; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + if (fmax_snan_l(q)) + { + q.m |= 1ULL << 62; + __builtin_memcpy(&y, &q, sizeof q); + return y; + } + if (fmax_nan_l(p)) + { + return fmax_nan_l(q) ? x : y; + } + if (fmax_nan_l(q)) + { + return x; + } + + if (x < y) + { + return y; + } + if (y < x) + { + return x; + } + if (x == 0.0L) + { + xneg = (int)(p.se >> 15); + yneg = (int)(q.se >> 15); + if (xneg != yneg) + { + return xneg ? y : x; + } + } + return x; +} + +/* + * The larger of two floats. + */ +float +fmaxf(float x, float y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmax_mm_f(x, y); +} + +/* + * The larger of two doubles. + */ +double +fmax(double x, double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmax_mm_d(x, y); +} + +/* + * The larger of two long doubles. + */ +long double +fmaxl(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmax_mm_l(x, y); +} diff --git a/src/math/fmin.c b/src/math/fmin.c new file mode 100644 index 0000000..69a3dd6 --- /dev/null +++ b/src/math/fmin.c @@ -0,0 +1,267 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * The smaller of x and y (C23 7.12.12.4), all three precisions. + * + * A quiet NaN argument is ignored in favor of the numeric one; a signaling + * NaN argument makes the result a quiet NaN (the host glibc returns the + * signaling argument quieted, payload preserved). When the arguments are + * +0.0 and -0.0 the result is -0.0, and a pair of equal nonzero values + * returns either. Numeric ordering is an ordinary comparison on the + * already-NaN-screened operands, so no exception is raised and the chosen + * operand is returned with its bits untouched. The NaN detection and the + * quiet splice are done on the raw IEEE 754 word (the __builtin_isnan + * classification would be equally exact but the word form also exposes the + * signaling bit); no rounding mode and no arithmetic instruction is + * involved, so the function is pure and fold-free at every level. + */ + +/* + * The double format: bit 63 is the sign, bits 62..52 the exponent biased + * by 1023, bits 51..0 the fraction; exponent 0x7ff with a nonzero fraction + * is a NaN, and the quiet bit is 0x0008000000000000 (bit 51). + */ +static int +fmin_nan_d(unsigned long long b) +{ + return ((b >> 52) & 0x7ff) == 0x7ff && (b & 0xFFFFFFFFFFFFFULL) != 0; +} + +static int +fmin_snan_d(unsigned long long b) +{ + return fmin_nan_d(b) && (b & 0x0008000000000000ULL) == 0; +} + +static double +fmin_mm_d(double x, double y, int want_max) +{ + unsigned long long xb; + unsigned long long yb; + int xneg; + int yneg; + + __builtin_memcpy(&xb, &x, sizeof xb); + __builtin_memcpy(&yb, &y, sizeof yb); + + /* A signaling NaN in either argument wins, quieted with its payload. */ + if (fmin_snan_d(xb)) + { + xb |= 0x0008000000000000ULL; + __builtin_memcpy(&x, &xb, sizeof xb); + return x; + } + if (fmin_snan_d(yb)) + { + yb |= 0x0008000000000000ULL; + __builtin_memcpy(&y, &yb, sizeof yb); + return y; + } + + /* Quiet NaNs are ignored: the numeric argument (or x, for two) wins. */ + if (fmin_nan_d(xb)) + { + return fmin_nan_d(yb) ? x : y; + } + if (fmin_nan_d(yb)) + { + return x; + } + + if (x < y) + { + return x; + } + if (y < x) + { + return y; + } + + /* Equal, so possibly a +0/-0 pair: fmin prefers -0, fmax prefers +0. */ + if (x == 0.0) + { + xneg = (int)(xb >> 63); + yneg = (int)(yb >> 63); + if (xneg != yneg) + { + return want_max ? (xneg ? y : x) : (xneg ? x : y); + } + } + return x; +} + +/* + * The float format: bit 31 the sign, bits 30..23 the exponent biased by + * 127, bits 22..0 the fraction; quiet bit 0x00400000 (bit 22). + */ +static int +fmin_nan_f(unsigned int b) +{ + return ((b >> 23) & 0xff) == 0xff && (b & 0x7FFFFFU) != 0; +} + +static int +fmin_snan_f(unsigned int b) +{ + return fmin_nan_f(b) && (b & 0x00400000U) == 0; +} + +static float +fmin_mm_f(float x, float y, int want_max) +{ + unsigned int xb; + unsigned int yb; + int xneg; + int yneg; + + __builtin_memcpy(&xb, &x, sizeof xb); + __builtin_memcpy(&yb, &y, sizeof yb); + + if (fmin_snan_f(xb)) + { + xb |= 0x00400000U; + __builtin_memcpy(&x, &xb, sizeof xb); + return x; + } + if (fmin_snan_f(yb)) + { + yb |= 0x00400000U; + __builtin_memcpy(&y, &yb, sizeof yb); + return y; + } + if (fmin_nan_f(xb)) + { + return fmin_nan_f(yb) ? x : y; + } + if (fmin_nan_f(yb)) + { + return x; + } + + if (x < y) + { + return x; + } + if (y < x) + { + return y; + } + if (x == 0.0f) + { + xneg = (int)(xb >> 31); + yneg = (int)(yb >> 31); + if (xneg != yneg) + { + return want_max ? (xneg ? y : x) : (xneg ? x : y); + } + } + return x; +} + +/* + * The x86 80-bit extended format: 64-bit significand m in bytes 0..7 (the + * integer bit is explicit) and a sign/exponent word se in bytes 8..9, sign + * in bit 15 and exponent biased by 16383 in bits 14..0. A NaN has + * (se & 0x7fff) == 0x7fff with a nonzero fraction; the quiet bit is bit 62 + * of m. + */ +struct fmin_ld_word +{ + unsigned long long m; + unsigned short se; +}; + +static int +fmin_nan_l(struct fmin_ld_word p) +{ + return (p.se & 0x7fff) == 0x7fff && (p.m & 0x7FFFFFFFFFFFFFFFULL) != 0; +} + +static int +fmin_snan_l(struct fmin_ld_word p) +{ + return fmin_nan_l(p) && (p.m & (1ULL << 62)) == 0; +} + +static long double +fmin_mm_l(long double x, long double y, int want_max) +{ + struct fmin_ld_word p; + struct fmin_ld_word q; + int xneg; + int yneg; + + __builtin_memcpy(&p, &x, sizeof p); + __builtin_memcpy(&q, &y, sizeof q); + + if (fmin_snan_l(p)) + { + p.m |= 1ULL << 62; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + if (fmin_snan_l(q)) + { + q.m |= 1ULL << 62; + __builtin_memcpy(&y, &q, sizeof q); + return y; + } + if (fmin_nan_l(p)) + { + return fmin_nan_l(q) ? x : y; + } + if (fmin_nan_l(q)) + { + return x; + } + + if (x < y) + { + return x; + } + if (y < x) + { + return y; + } + if (x == 0.0L) + { + xneg = (int)(p.se >> 15); + yneg = (int)(q.se >> 15); + if (xneg != yneg) + { + return want_max ? (xneg ? y : x) : (xneg ? x : y); + } + } + return x; +} + +/* + * The smaller of two floats. + */ +float +fminf(float x, float y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmin_mm_f(x, y, 0); +} + +/* + * The smaller of two doubles. + */ +double +fmin(double x, double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmin_mm_d(x, y, 0); +} + +/* + * The smaller of two long doubles. + */ +long double +fminl(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmin_mm_l(x, y, 0); +} diff --git a/src/math/fmod.c b/src/math/fmod.c new file mode 100644 index 0000000..9bdae86 --- /dev/null +++ b/src/math/fmod.c @@ -0,0 +1,423 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#ifdef HAVE_CONFIG_H +#include +#endif + +/* + * The floating-point remainder x - n*y with n = x/y truncated toward zero + * (C23 7.12.10.1), all three precisions. The remainder is always exact + * and carries x's sign (a zero result is a signed zero with x's sign), so + * every implementation must agree bit for bit; the result is computed with + * a restoring long division on the significands, in the style of the musl + * fmod family but derived here from first principles. + * + * Each operand is normalized into + * + * value = m * 2^(ex - W) + * + * where m is the significand with its msb pinned to bit P (P = 52/23/63 + * and W = 1075/150/16446 for double/float/80-bit) and ex is the exponent + * in the same units as the biased field (so ex is the stored exponent for + * normals and goes negative for subnormals). Aligned this way, one binary + * long-division step per exponent difference decides whether a multiple of + * the divisor fits: the current remainder significand is compared against + * the divisor significand, the divisor is subtracted once when it fits, + * and the remainder is doubled for the next, half-weight, digit. The + * comparison invariant keeps the remainder below twice the divisor at + * every step, so one subtraction per digit is always enough and every + * subtraction (and the exact zero test) is a plain integer operation. + * + * Domain errors mirror the host glibc: fmod(+-0, +-0), fmod(+-Inf, y) and + * any NaN argument return a NaN (the classic indefinite pattern, which is + * what glibc's x87 path produces) and, in the library build only, set + * errno to EDOM. fmod(x, +-Inf) is x for a finite x, and |x| <= |y| hands + * x back unchanged, which keeps the signed-zero and exact cases exact. + */ + +/* The double format: bit 63 the sign, bits 62..52 the exponent biased by + * 1023, bits 51..0 the fraction. Value = m * 2^(ex - 1075) with the msb + * of m at bit 52. */ +static double +fmod_d(double x, double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + unsigned long long xb; + unsigned long long yb; + unsigned long long sx; + unsigned long long mx; + unsigned long long my; + unsigned long long s; + unsigned long long i; + int ex; + int ey; + + __builtin_memcpy(&xb, &x, sizeof xb); + __builtin_memcpy(&yb, &y, sizeof yb); + sx = xb & (1ULL << 63); + ex = (int)((xb >> 52) & 0x7ff); + ey = (int)((yb >> 52) & 0x7ff); + + /* x not finite, y zero, or y a NaN: domain error. */ + if (ex == 0x7ff || (yb << 1) == 0 || (ey == 0x7ff && (yb & 0xFFFFFFFFFFFFFULL) != 0)) + { +#ifdef HAVE_CONFIG_H + errno = EDOM; +#endif + return __builtin_nan(""); /* indefinite NaN, sign negative */ + } + /* x zero or a finite x below an infinite y: x unchanged. (x a NaN or + * +-Inf fell into the domain branch above, which returns a NaN for + * both.) */ + if ((xb << 1) == 0 || ey == 0x7ff) + { + return x; + } + + /* Compare magnitudes directly on the (biased exponent, fraction) + * words: finite IEEE magnitudes order lexicographically, so a plain + * integer compare of the sign-stripped words decides |x| vs |y|. */ + if ((xb & 0x7FFFFFFFFFFFFFFFULL) <= (yb & 0x7FFFFFFFFFFFFFFFULL)) + { + if ((xb & 0x7FFFFFFFFFFFFFFFULL) == (yb & 0x7FFFFFFFFFFFFFFFULL)) + { + return x * 0.0; /* exact: result is +-0 with x's sign */ + } + return x; + } + + /* Normalize x: subnormals (ex == 0) get their msb shifted up to bit + * 52 and ex counts the shift below the smallest normal. */ + mx = xb & 0xFFFFFFFFFFFFFULL; + if (ex == 0) + { + s = 63 - (unsigned long long)__builtin_clzll(mx); + mx <<= (52 - s); + ex = (int)s - 51; + } + else + { + mx |= 1ULL << 52; + } + my = yb & 0xFFFFFFFFFFFFFULL; + if (ey == 0) + { + s = 63 - (unsigned long long)__builtin_clzll(my); + my <<= (52 - s); + ey = (int)s - 51; + } + else + { + my |= 1ULL << 52; + } + + /* Long division: subtract the aligned divisor significand once per + * bit of quotient, doubling the remainder between bits. */ + for (; ex > ey; ex--) + { + if (mx >= my) + { + i = mx - my; + if (i == 0) + { + return x * 0.0; + } + mx = i; + } + mx <<= 1; + } + i = mx - my; + if (mx >= my) + { + if (i == 0) + { + return x * 0.0; + } + mx = i; + } + + /* Bring the msb of the remainder back to bit 52. */ + for (; (mx >> 52) == 0; mx <<= 1, ex--) + { + } + + if (ex > 0) + { + xb = sx | ((unsigned long long)ex << 52) | (mx - (1ULL << 52)); + } + else + { + xb = sx | (mx >> (1 - ex)); + } + __builtin_memcpy(&x, &xb, sizeof xb); + return x; +} + +/* The float format: bit 31 the sign, bits 30..23 the exponent biased by + * 127, bits 22..0 the fraction. Value = m * 2^(ex - 150) with the msb of + * m at bit 23. */ +static float +fmod_f(float x, float y) // NOLINT(bugprone-easily-swappable-parameters) +{ + unsigned int xb; + unsigned int yb; + unsigned int sx; + unsigned int mx; + unsigned int my; + unsigned int s; + unsigned int i; + int ex; + int ey; + + __builtin_memcpy(&xb, &x, sizeof xb); + __builtin_memcpy(&yb, &y, sizeof yb); + sx = xb & (1U << 31); + ex = (int)((xb >> 23) & 0xff); + ey = (int)((yb >> 23) & 0xff); + + if (ex == 0xff || (yb << 1) == 0 || (ey == 0xff && (yb & 0x7FFFFFU) != 0)) + { +#ifdef HAVE_CONFIG_H + errno = EDOM; +#endif + return __builtin_nanf(""); + } + if ((xb << 1) == 0 || ey == 0xff) + { + return x; + } + + if ((xb & 0x7FFFFFFFU) <= (yb & 0x7FFFFFFFU)) + { + if ((xb & 0x7FFFFFFFU) == (yb & 0x7FFFFFFFU)) + { + return x * 0.0f; + } + return x; + } + + mx = xb & 0x7FFFFFU; + if (ex == 0) + { + s = 31 - (unsigned int)__builtin_clz(mx); + mx <<= (23 - s); + ex = (int)s - 22; + } + else + { + mx |= 1U << 23; + } + my = yb & 0x7FFFFFU; + if (ey == 0) + { + s = 31 - (unsigned int)__builtin_clz(my); + my <<= (23 - s); + ey = (int)s - 22; + } + else + { + my |= 1U << 23; + } + + for (; ex > ey; ex--) + { + if (mx >= my) + { + i = mx - my; + if (i == 0) + { + return x * 0.0f; + } + mx = i; + } + mx <<= 1; + } + i = mx - my; + if (mx >= my) + { + if (i == 0) + { + return x * 0.0f; + } + mx = i; + } + + for (; (mx >> 23) == 0; mx <<= 1, ex--) + { + } + + if (ex > 0) + { + xb = sx | ((unsigned int)ex << 23) | (mx - (1U << 23)); + } + else + { + xb = sx | (mx >> (1 - ex)); + } + __builtin_memcpy(&x, &xb, sizeof xb); + return x; +} + +/* + * The x86 80-bit extended format: 64-bit significand m with an explicit + * integer bit in bytes 0..7, sign/exponent word se in bytes 8..9. Value = + * m * 2^(ex - 16446) with the msb of m at bit 63; a canonical normal has + * ex = se & 0x7fff, and a subnormal (se field 0) or a defensive unnormal + * (nonzero field with m < 2^63) is normalized by shifting m's msb up to + * bit 63. The division reuses the same digit loop as the narrower + * formats; only the doubling step differs, because a full-width + * significand has no free bit above its msb (see the three-way branch in + * the loop below). + */ +struct fmod_ld_word +{ + unsigned long long m; + unsigned short se; +}; + +static long double +fmod_l(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + struct fmod_ld_word p; + struct fmod_ld_word q; + unsigned long long s; + unsigned long long mx; + unsigned long long my; + unsigned long long i; + int ex; + int ey; + + __builtin_memcpy(&p, &x, sizeof p); + __builtin_memcpy(&q, &y, sizeof q); + ex = p.se & 0x7fff; + ey = q.se & 0x7fff; + + /* x not finite, y zero, or y a NaN: domain error. */ + if (ex == 0x7fff || (q.m == 0 && ey == 0) || + (ey == 0x7fff && (q.m & 0x7FFFFFFFFFFFFFFFULL) != 0)) + { +#ifdef HAVE_CONFIG_H + errno = EDOM; +#endif + p.m = 0xC000000000000000ULL; + p.se = 0xFFFF; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + /* x zero or a finite x below an infinite y: x unchanged. (x a NaN or + * +-Inf fell into the domain branch above.) */ + if ((p.m == 0 && ex == 0) || ey == 0x7fff) + { + return x; + } + + /* Normalize both significands to the [2^63, 2^64) frame; a subnormal + * (field 0) or defensive unnormal (field nonzero, msb below 63) has + * its msb shifted up to bit 63 with the exponent adjusted, so every + * value below obeys value = m * 2^(ex - 16446). */ + if (p.m < 0x8000000000000000ULL) + { + s = 63 - (unsigned long long)__builtin_clzll(p.m); + p.m <<= (63 - s); + ex = (ex == 0 ? 1 : ex) + (int)s - 63; + } + if (q.m < 0x8000000000000000ULL) + { + s = 63 - (unsigned long long)__builtin_clzll(q.m); + q.m <<= (63 - s); + ey = (ey == 0 ? 1 : ey) + (int)s - 63; + } + mx = p.m; + my = q.m; + + if (ex < ey || (ex == ey && mx <= my)) + { + if (ex == ey && mx == my) + { + return x * 0.0L; + } + return x; + } + + /* Long division with a full-width significand: a 64-bit divisor + * leaves no headroom above its own msb for the per-bit doubling, so + * the doubled remainder that overflows is exactly one divisor at the + * next, half, weight and is absorbed by a subtraction there. */ + for (; ex > ey; ex--) + { + i = mx - my; + if (mx >= my) + { + if (i == 0) + { + return x * 0.0L; + } + mx = 2 * i; + } + else if (2 * mx < mx) + { + mx = 2 * mx - my; + } + else + { + mx = 2 * mx; + } + } + i = mx - my; + if (mx >= my) + { + if (i == 0) + { + return x * 0.0L; + } + mx = i; + } + + /* Bring the msb of the remainder back to bit 63. */ + for (; mx < 0x8000000000000000ULL; mx <<= 1, ex--) + { + } + + if (ex > 0) + { + p.se = (unsigned short)((p.se & 0x8000) | (unsigned short)ex); + p.m = mx; + } + else + { + p.se = (unsigned short)(p.se & 0x8000); + p.m = mx >> (1 - ex); + } + __builtin_memcpy(&x, &p, sizeof p); + return x; +} + +/* + * The float remainder, as fmod_d. + */ +float +fmodf(float x, float y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmod_f(x, y); +} + +/* + * The double remainder, as fmod_d. + */ +double +fmod(double x, double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmod_d(x, y); +} + +/* + * The long double remainder, as fmod_d. + */ +long double +fmodl(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return fmod_l(x, y); +} diff --git a/src/math/frexp.c b/src/math/frexp.c new file mode 100644 index 0000000..703e475 --- /dev/null +++ b/src/math/frexp.c @@ -0,0 +1,173 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Split x into a fraction f in [1/2, 1) (or 0) and an integer exponent + * *exp with x == f * 2^(*exp), all three precisions (C23 7.12.6.4). + * frexp(±0) returns x with *exp 0; frexp(±Inf) and frexp(NaN) return x + * with *exp 0 (the glibc behavior the tests pin down; the C standard + * leaves the Inf/NaN exponent unspecified). The fraction keeps x's sign, + * so frexp(-6.0) is -0.75 with *exp 3 and frexp(0x1p-1074) is 0.5 with + * *exp -1073. + * + * Each precision works from the same unified view used across the + * src/math/ slices (see math_impl.h): a finite value is m * 2^p with m the + * integer significand carrying the explicit integer bit. The fraction is + * m * 2^-w (w = 53/24/64 significand bits), which sits in [1/2, 1), so + * *exp is p + w and the result keeps the significand's low bits — for a + * normal input the fraction field is simply unchanged and only the + * exponent field is rewritten. + */ + +/* + * The double format: bit 63 the sign, bits 62..52 the exponent biased by + * 1023, bits 51..0 the fraction, 53-bit significand. A normal input has + * p = ef - 1075, so *exp = ef - 1022 and the fraction field is the input's + * own 52-bit fraction (the implicit bit stays put as the leading 1 of the + * fraction's [1/2, 1) significand). A subnormal input with fraction msb + * at position s (0..51) is value frac * 2^-1074; shifting frac up to the + * implicit-bit position makes m = frac << (52 - s), hence *exp = s - 1073 + * and fraction field m - 2^52. + */ +double +frexp(double x, int *exp) +{ + unsigned long long bits; + unsigned long long sign; + unsigned long long ef; + unsigned long long frac; + unsigned long long m; + int s; + + __builtin_memcpy(&bits, &x, sizeof bits); + sign = bits & (1ULL << 63); + ef = (bits >> 52) & 0x7ff; + frac = bits & 0xFFFFFFFFFFFFFULL; + + if (ef == 0x7ff) + { + *exp = 0; + return x; /* ±Inf and NaN */ + } + if (ef == 0) + { + if (frac == 0) + { + *exp = 0; + return x; /* ±0 */ + } + s = 63 - __builtin_clzll(frac); + m = frac << (52 - s); + *exp = s - 1073; + bits = sign | (0x3feULL << 52) | (m - (1ULL << 52)); + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + + *exp = (int)ef - 1022; + bits = sign | (0x3feULL << 52) | frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x; +} + +/* + * The float format: bit 31 the sign, bits 30..23 the exponent biased by + * 127, bits 22..0 the fraction, 24-bit significand. Normal *exp = + * ef - 126; a subnormal input with fraction msb at s (0..22) normalizes to + * m = frac << (23 - s) and has *exp = s - 148. + */ +float +frexpf(float x, int *exp) +{ + unsigned int bits; + unsigned int sign; + unsigned int ef; + unsigned int frac; + unsigned int m; + int s; + + __builtin_memcpy(&bits, &x, sizeof bits); + sign = bits & (1U << 31); + ef = (bits >> 23) & 0xff; + frac = bits & 0x7FFFFFU; + + if (ef == 0xff) + { + *exp = 0; + return x; + } + if (ef == 0) + { + if (frac == 0) + { + *exp = 0; + return x; + } + s = 31 - __builtin_clz(frac); + m = frac << (23 - s); + *exp = s - 148; + bits = sign | (0x7eU << 23) | (m - (1U << 23)); + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + + *exp = (int)ef - 126; + bits = sign | (0x7eU << 23) | frac; + __builtin_memcpy(&x, &bits, sizeof x); + return x; +} + +/* + * The x86 80-bit extended format: 64 significand bits m with an explicit + * integer bit, plus a sign/exponent word se biased by 16383. A canonical + * normal has m in [2^63, 2^64) and p = ef - 16446, so *exp = ef - 16382 + * and the fraction is m * 2^-64 with se rewritten to 16382. A subnormal + * (ef == 0) is value m * 2^-16445; shifting m's msb (position s, 0..62) + * up to bit 63 gives *exp = s - 16444. Unnormal inputs (ef > 0 with + * m < 2^63) normalize the same way with *exp = ef + s - 16445. + */ +long double +frexpl(long double x, int *exp) +{ + struct + { + unsigned long long m; + unsigned short se; + } p; + unsigned long long m; + int ef; + int s; + + __builtin_memcpy(&p, &x, sizeof p); + m = p.m; + ef = p.se & 0x7fff; + + if (ef == 0x7fff) + { + *exp = 0; + return x; + } + if (m == 0) + { + *exp = 0; + return x; /* ±0 and empty degenerate encodings */ + } + if (m < 0x8000000000000000ULL) + { + s = 63 - __builtin_clzll(m); + m <<= (63 - s); + *exp = (ef == 0 ? 1 : ef) + s - 16445; + p.se = (unsigned short)((p.se & 0x8000) | 0x3ffe); + p.m = m; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + + *exp = ef - 16382; + p.se = (unsigned short)((p.se & 0x8000) | 0x3ffe); + __builtin_memcpy(&x, &p, sizeof p); + return x; +} diff --git a/src/math/ilogb.c b/src/math/ilogb.c new file mode 100644 index 0000000..a07f92e --- /dev/null +++ b/src/math/ilogb.c @@ -0,0 +1,126 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * The signed exponent of x as an int (C23 7.12.6.5p1), all three + * precisions: ilogb(x) is floor(log2 |x|) for a nonzero finite x. The + * sentinel returns are FP_ILOGB0 (INT_MIN) for ±0 and FP_ILOGBNAN + * (INT_MAX) for ±Inf and NaN, without touching errno, so the functions + * carry the const attribute and never trap. Subnormal arguments are + * normalized by scanning the fraction for its highest set bit. + * + * ilogb(1.0) == 0, ilogb(8.0) == 3, ilogb(0.5) == -1, and the smallest + * subnormal double, 2^-1074, yields -1074. + */ + +/* + * The double format: bits 62..52 hold the exponent biased by 1023 and + * bits 51..0 the fraction; a subnormal with fraction F (no implicit + * bit) is F * 2^-1074, so its exponent is the fraction's highest set + * bit position minus 1074. + */ +int +ilogb(double x) +{ + unsigned long long bits; + unsigned long long frac; + int e; + int msb; + + __builtin_memcpy(&bits, &x, sizeof bits); + frac = bits & 0xFFFFFFFFFFFFFULL; + e = (int)((bits >> 52) & 0x7ff); + + if (e == 0x7ff) + { + /* ±Inf and NaN both map to FP_ILOGBNAN. */ + return FP_ILOGBNAN; + } + if (e == 0) + { + if (frac == 0) + { + return FP_ILOGB0; + } + msb = 63 - __builtin_clzll(frac); + return msb - 1074; + } + return e - 1023; +} + +/* + * The float format: bits 30..23 hold the exponent biased by 127 and + * bits 22..0 the fraction; a subnormal with fraction F is F * 2^-149. + */ +int +ilogbf(float x) +{ + unsigned int bits; + unsigned int frac; + int e; + int msb; + + __builtin_memcpy(&bits, &x, sizeof bits); + frac = bits & 0x7FFFFFU; + e = (int)((bits >> 23) & 0xff); + + if (e == 0xff) + { + return FP_ILOGBNAN; + } + if (e == 0) + { + if (frac == 0) + { + return FP_ILOGB0; + } + msb = 31 - __builtin_clz(frac); + return msb - 149; + } + return e - 127; +} + +/* + * The x86 80-bit extended format: 64 significand bits m (explicit + * integer bit) with a sign/exponent word se biased by 16383. A normal m + * in [2^63, 2^64) sits in [2^(e-16383), 2^(e-16382)). A subnormal with + * e == 0 is interpreted as if e were 1, i.e. m * 2^-16445; an unnormal + * (e > 0 with m < 2^63) is m * 2^(e - 16446). + */ +int +ilogbl(long double x) +{ + struct + { + unsigned long long m; + unsigned short se; + } p; + int e; + int msb; + + __builtin_memcpy(&p, &x, sizeof p); + e = p.se & 0x7fff; + + if (e == 0x7fff) + { + return FP_ILOGBNAN; + } + if (p.m == 0) + { + return FP_ILOGB0; + } + if (e == 0) + { + msb = 63 - __builtin_clzll(p.m); + return msb - 16445; + } + if (p.m < 0x8000000000000000ULL) + { + msb = 63 - __builtin_clzll(p.m); + return msb + e - 16446; + } + return e - 16383; +} diff --git a/src/math/ldexp.c b/src/math/ldexp.c new file mode 100644 index 0000000..7a1825a --- /dev/null +++ b/src/math/ldexp.c @@ -0,0 +1,76 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "math_impl.h" + +/* + * x times 2^n (C23 7.12.6.3), all three precisions: the inverse of frexp. + * Scaling a finite value by a power of two is exact whenever the result is + * a normal number; only a result that lands in the subnormal range can + * need rounding (the vl_rne_u64 round-to-nearest-even shift in the shared + * scale cores, see math_impl.h). Overflow returns +-Inf and, in the + * library build, sets errno to ERANGE as POSIX requires. + * + * The errno write is a TCB dereference that only exists when this file is + * compiled as part of the real library (config.h present). The host- + * linked standalone test binaries compile these sources without + * HAVE_CONFIG_H, so they never reference __errno_location and never touch + * the host's errno thread slot; the tests therefore check values only. + * ldexp(±0) is ±0 for any n, and ±Inf/NaN pass through unchanged. + */ + +/* + * Scale a double: the shared core reports overflow through the flag and + * this wrapper raises ERANGE when the library build demands it. + */ +double +ldexp(double x, int n) +{ + int overflowed = 0; + double r = vl_scale2_d(x, n, &overflowed); + +#ifdef HAVE_CONFIG_H + if (overflowed) + { + errno = ERANGE; + } +#endif + return r; +} + +float +ldexpf(float x, int n) +{ + int overflowed = 0; + float r = vl_scale2_f(x, n, &overflowed); + +#ifdef HAVE_CONFIG_H + if (overflowed) + { + errno = ERANGE; + } +#endif + return r; +} + +long double +ldexpl(long double x, int n) +{ + int overflowed = 0; + long double r = vl_scale2_ld(x, n, &overflowed); + +#ifdef HAVE_CONFIG_H + if (overflowed) + { + errno = ERANGE; + } +#endif + return r; +} diff --git a/src/math/llrint.c b/src/math/llrint.c new file mode 100644 index 0000000..e015c41 --- /dev/null +++ b/src/math/llrint.c @@ -0,0 +1,45 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * The nearest integral value to x in the current rounding direction, + * returned as long long (C23 7.12.9.8), all three precisions. A result + * outside the range of long long is a range error whose return value is + * unspecified, so the header declares no const attribute; the test + * corpus keeps |x| < 2^62 where every result is exact. + * + * As for lrint (see lrint.c): the __builtin_llrint* forms are never + * folded on this target, so llrint is built as rint-then-convert, with + * __builtin_rint* folding to the in-line round-to-nearest-even sequence + * and the cast to long long being exact on the integral result. + */ + +/* + * As llrint, for a float argument. + */ +long long +llrintf(float x) +{ + return (long long)__builtin_rintf(x); +} + +/* + * As llrint, for a double argument. + */ +long long +llrint(double x) +{ + return (long long)__builtin_rint(x); +} + +/* + * As llrint, for a long double argument. + */ +long long +llrintl(long double x) +{ + return (long long)__builtin_rintl(x); +} diff --git a/src/math/llround.c b/src/math/llround.c new file mode 100644 index 0000000..17725b5 --- /dev/null +++ b/src/math/llround.c @@ -0,0 +1,73 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +/* + * The nearest integral value to x, halfway cases rounded away from zero, + * returned as long long (C23 7.12.9.10), all three precisions. + * + * Rounding is delegated to the round/roundf/roundl implementations (see + * round.c): on this target every __builtin_round* is an external call even + * at -O2 -- no SSE4.1 in the default -march -- so round() here reaches the + * library's own half-away-from-zero rounding, returning an exactly + * integral value of the same floating type. Narrowing that value to long + * long is then a plain, exact conversion whenever it is in range. + * + * The unrepresentable cases mirror the host glibc on x86-64, whose lround + * family is the hardware cvttsd2si sequence: a NaN, an Inf, and any finite + * result of magnitude >= 2^63 (the smallest such rounded value is exactly + * 2^63) all collapse to LLONG_MIN with errno untouched -- measured against + * glibc 2.44, which returns LLONG_MIN for 1e300, -1e300, +-Inf and NaN and + * never sets errno. The range checks happen in the argument's own + * precision (0x1p63 is exact in float, double and the x87 extended + * format), and -2^63 -- exactly representable and equal to LLONG_MIN -- is + * deliberately allowed through the strict lower-bound test. + */ + +/* + * As llround, for a float argument. + */ +long long +llroundf(float x) +{ + float r = roundf(x); + + if (isnan(r) || r >= 0x1p63f || r < -0x1p63f) + { + return LLONG_MIN; + } + return (long long)r; +} + +/* + * As llround, for a double argument. + */ +long long +llround(double x) +{ + double r = round(x); + + if (isnan(r) || r >= 0x1p63 || r < -0x1p63) + { + return LLONG_MIN; + } + return (long long)r; +} + +/* + * As llround, for a long double argument. + */ +long long +llroundl(long double x) +{ + long double r = roundl(x); + + if (isnan(r) || r >= 0x1p63L || r < -0x1p63L) + { + return LLONG_MIN; + } + return (long long)r; +} diff --git a/src/math/logb.c b/src/math/logb.c new file mode 100644 index 0000000..2256984 --- /dev/null +++ b/src/math/logb.c @@ -0,0 +1,123 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * The signed exponent of x as a floating-point value (C23 7.12.6.6), + * all three precisions: logb(x) is floor(log2 |x|) as a float/double/ + * long double. logb(±0) is -Inf and logb(±Inf) is +Inf, each raising + * the division-by-zero/invalid exception through the hardware but + * without an errno path in the representable domain; logb(NaN) returns + * the NaN. Pure, so the header marks the functions const; the exponent + * extraction is exactly the ilogb one (see ilogb.c) with the int result + * converted back to the argument's precision. + */ + +/* + * The double format: bits 62..52 hold the exponent biased by 1023 and + * bits 51..0 the fraction. A subnormal with fraction F is F * 2^-1074. + */ +double +logb(double x) +{ + unsigned long long bits; + unsigned long long frac; + int e; + int msb; + + __builtin_memcpy(&bits, &x, sizeof bits); + frac = bits & 0xFFFFFFFFFFFFFULL; + e = (int)((bits >> 52) & 0x7ff); + + if (e == 0x7ff) + { + /* Inf maps to +Inf (both signs); NaN passes through. */ + return frac == 0 ? HUGE_VAL : x; + } + if (e == 0) + { + if (frac == 0) + { + return -HUGE_VAL; + } + msb = 63 - __builtin_clzll(frac); + return (double)(msb - 1074); + } + return (double)(e - 1023); +} + +/* + * The float format: bits 30..23 hold the exponent biased by 127 and + * bits 22..0 the fraction; a subnormal with fraction F is F * 2^-149. + */ +float +logbf(float x) +{ + unsigned int bits; + unsigned int frac; + int e; + int msb; + + __builtin_memcpy(&bits, &x, sizeof bits); + frac = bits & 0x7FFFFFU; + e = (int)((bits >> 23) & 0xff); + + if (e == 0xff) + { + return frac == 0 ? HUGE_VALF : x; + } + if (e == 0) + { + if (frac == 0) + { + return -HUGE_VALF; + } + msb = 31 - __builtin_clz(frac); + return (float)(msb - 149); + } + return (float)(e - 127); +} + +/* + * The x86 80-bit extended format: 64 significand bits m (explicit + * integer bit) with a sign/exponent word se biased by 16383. A normal m + * in [2^63, 2^64) sits in [2^(e-16383), 2^(e-16382)). A subnormal with + * e == 0 is interpreted as if e were 1, i.e. m * 2^-16445; an unnormal + * (e > 0 with m < 2^63) is m * 2^(e - 16446). + */ +long double +logbl(long double x) +{ + struct + { + unsigned long long m; + unsigned short se; + } p; + int e; + int msb; + + __builtin_memcpy(&p, &x, sizeof p); + e = p.se & 0x7fff; + + if (e == 0x7fff) + { + return (p.m & 0x7FFFFFFFFFFFFFFFULL) == 0 ? HUGE_VALL : x; + } + if (p.m == 0) + { + return -HUGE_VALL; + } + if (e == 0) + { + msb = 63 - __builtin_clzll(p.m); + return (long double)(msb - 16445); + } + if (p.m < 0x8000000000000000ULL) + { + msb = 63 - __builtin_clzll(p.m); + return (long double)(msb + e - 16446); + } + return (long double)(e - 16383); +} diff --git a/src/math/lrint.c b/src/math/lrint.c new file mode 100644 index 0000000..736e3c1 --- /dev/null +++ b/src/math/lrint.c @@ -0,0 +1,49 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * The nearest integral value to x in the current rounding direction, + * returned as long (C23 7.12.9.7), all three precisions. A result + * outside the range of long is a range error whose return value is + * unspecified, so the header declares no const attribute; the test + * corpus keeps |x| < 2^62 where every result is exact. + * + * GCC never folds the __builtin_lrint* forms on this target (external + * lrint@PLT calls at every optimization level), so lrint is built as + * rint-then-convert: __builtin_rint* folds to the in-line round-to- + * nearest-even sequence (see rint.c), producing an exact integral value, + * and the cast to long is then exact no matter which conversion + * instruction GCC emits. The rint step honors the MXCSR/x87 rounding + * mode, the only reachable one being the default round-to-nearest-even + * (no exists in vlibc yet). + */ + +/* + * As lrint, for a float argument. + */ +long +lrintf(float x) +{ + return (long)__builtin_rintf(x); +} + +/* + * As lrint, for a double argument. + */ +long +lrint(double x) +{ + return (long)__builtin_rint(x); +} + +/* + * As lrint, for a long double argument. + */ +long +lrintl(long double x) +{ + return (long)__builtin_rintl(x); +} diff --git a/src/math/lround.c b/src/math/lround.c new file mode 100644 index 0000000..087ac98 --- /dev/null +++ b/src/math/lround.c @@ -0,0 +1,64 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +/* + * The nearest integral value to x, halfway cases rounded away from zero, + * returned as long (C23 7.12.9.9), all three precisions. + * + * On this LP64 target long is 64 bits wide, so lround is llround in every + * observable way: the value is rounded half away from zero by round.c and + * narrowed, and a NaN, an Inf, or any rounded magnitude >= 2^63 collapses + * to LONG_MIN, exactly what the host glibc's hardware-conversion lround + * returns for such arguments (measured on glibc 2.44; errno untouched). + * The reasoning, the range test shape and the -2^63 == LONG_MIN boundary + * subtlety are those documented in llround.c. + */ + +/* + * As lround, for a float argument. + */ +long +lroundf(float x) +{ + float r = roundf(x); + + if (isnan(r) || r >= 0x1p63f || r < -0x1p63f) + { + return LONG_MIN; + } + return (long)r; +} + +/* + * As lround, for a double argument. + */ +long +lround(double x) +{ + double r = round(x); + + if (isnan(r) || r >= 0x1p63 || r < -0x1p63) + { + return LONG_MIN; + } + return (long)r; +} + +/* + * As lround, for a long double argument. + */ +long +lroundl(long double x) +{ + long double r = roundl(x); + + if (isnan(r) || r >= 0x1p63L || r < -0x1p63L) + { + return LONG_MIN; + } + return (long)r; +} diff --git a/src/math/math_impl.h b/src/math/math_impl.h new file mode 100644 index 0000000..1a04e2b --- /dev/null +++ b/src/math/math_impl.h @@ -0,0 +1,341 @@ +#ifndef VLIBC_MATH_IMPL_H +#define VLIBC_MATH_IMPL_H + +#include + +/* + * vlibc — private helpers shared by the src/math/ implementation files. + * + * This header is internal to the todo-39 arithmetic slices (frexp, ldexp, + * modf, scalbn, scalbln, and later slices that need to scale a value by a + * power of two); it is never installed and is not public API. + * + * Every real function below is spelled for all three precisions around one + * unified representation of a finite value: + * + * value = m * 2^p + * + * where m is an integer significand that carries the explicit integer bit + * (m in [2^52, 2^53) for double, [2^23, 2^24) for float, [2^63, 2^64) for + * the x86 80-bit extended format) and p is the exact unbiased power of + * two. Normalizing a subnormal input (and an 80-bit unnormal) into this + * shape is an exact left shift, so the scaling logic that follows never + * has to special-case the input class again. + * + * Scaling by 2^n only moves p: value = m * 2^(p + n). The result class + * is read straight off the target exponent k = p + n, and a right shift + * with round-to-nearest-even (vl_rne_u64) is needed only when the scaled + * value lands in the subnormal range and some low significand bits must be + * dropped. + * + * The format constants (mask widths, exponent biases) are the ones the + * rounding-family files (round.c, rint.c, ilogb.c, ...) already document; + * only the value interpretation changes here. + */ + +/* + * Round m >> r to nearest, ties to even, returning the kept integer. r is + * in [1, 64]; for r == 64 the entire significand is dropped and the only + * values that can survive are those above the tie point (1 is returned), + * while r > 64 can never round up because m < 2^64 is below the half-way + * threshold. Callers guarantee m < 2^64 and r >= 1. + */ +static inline unsigned long long +vl_rne_u64(unsigned long long m, int r) +{ + unsigned long long kept; + unsigned long long dropped; + unsigned long long half; + + if (r >= 64) + { + return (r > 64 || m <= (1ULL << 63)) ? 0 : 1; + } + kept = m >> r; + dropped = m & ((1ULL << r) - 1ULL); + half = 1ULL << (r - 1); + if (dropped > half || (dropped == half && (kept & 1ULL) != 0)) + { + kept++; + } + return kept; +} + +/* + * Scale a double by 2^n. Bits 62..52 are the exponent biased by 1023, + * bits 51..0 the fraction; the significand is 53 bits wide, so a normal + * input has m = 2^52 | frac and p = ef - 1075. A subnormal (ef == 0, + * frac != 0) with msb at position s (0..51) normalizes exactly to + * m = frac << (52 - s) with p = s - 1126. + * + * Target k = p + n: k >= 972 overflows (the exact result exceeds DBL_MAX; + * the flag is raised and +-Inf returned), k in [-1074, 971] is an exact + * normal result, and k <= -1075 is subnormal-or-zero with r = -1074 - k + * dropped bits; r >= 54 always rounds to zero, and a rounded-up field of + * exactly 2^52 is the smallest normal (ef 1, frac 0), reached when the + * exact value is the tie just above the largest subnormal. + */ +static inline double +vl_scale2_d(double x, int n, int *overflowed) // NOLINT(bugprone-easily-swappable-parameters) +{ + unsigned long long bits; + unsigned long long sign; + unsigned long long ef; + unsigned long long frac; + unsigned long long m; + unsigned long long f; + int s; + int p; + int k; + int r; + + __builtin_memcpy(&bits, &x, sizeof bits); + sign = bits & (1ULL << 63); + ef = (bits >> 52) & 0x7ff; + frac = bits & 0xFFFFFFFFFFFFFULL; + + if (ef == 0x7ff) + { + return x; /* +-Inf and NaN pass through unchanged */ + } + if (ef == 0) + { + if (frac == 0) + { + return x; /* +-0 stays +-0 for any n */ + } + s = 63 - __builtin_clzll(frac); + m = frac << (52 - s); + p = s - 1126; + } + else + { + m = frac | (1ULL << 52); + p = (int)ef - 1075; + } + + k = p + n; + if (k >= 972) + { + if (overflowed != NULL) + { + *overflowed = 1; + } + bits = sign | (0x7ffULL << 52); + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + if (k >= -1074) + { + bits = sign | ((unsigned long long)(k + 1075) << 52) | (m - (1ULL << 52)); + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + + r = -1074 - k; + if (r >= 54) + { + bits = sign; + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + f = vl_rne_u64(m, r); + if (f == 0) + { + bits = sign; + } + else if (f == (1ULL << 52)) + { + bits = sign | (1ULL << 52); /* smallest normal, reached by rounding up */ + } + else + { + bits = sign | f; + } + __builtin_memcpy(&x, &bits, sizeof x); + return x; +} + +/* + * Scale a float by 2^n. Bits 30..23 are the exponent biased by 127, bits + * 22..0 the fraction; a normal input has m = 2^23 | frac and p = ef - 150. + * A subnormal (ef == 0, frac != 0) with msb at s (0..22) normalizes to + * m = frac << (23 - s) with p = s - 172. + * + * k >= 105 overflows, k in [-149, 104] is an exact normal result, and + * k <= -150 is subnormal-or-zero with r = -149 - k; r >= 25 always rounds + * to zero and a rounded-up field of exactly 2^23 is the smallest normal. + */ +static inline float +vl_scale2_f(float x, int n, int *overflowed) // NOLINT(bugprone-easily-swappable-parameters) +{ + unsigned int bits; + unsigned int sign; + unsigned int ef; + unsigned int frac; + unsigned int m; + unsigned int f; + int s; + int p; + int k; + int r; + + __builtin_memcpy(&bits, &x, sizeof bits); + sign = bits & (1U << 31); + ef = (bits >> 23) & 0xff; + frac = bits & 0x7FFFFFU; + + if (ef == 0xff) + { + return x; + } + if (ef == 0) + { + if (frac == 0) + { + return x; + } + s = 31 - __builtin_clz(frac); + m = frac << (23 - s); + p = s - 172; + } + else + { + m = frac | (1U << 23); + p = (int)ef - 150; + } + + k = p + n; + if (k >= 105) + { + if (overflowed != NULL) + { + *overflowed = 1; + } + bits = sign | (0xffU << 23); + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + if (k >= -149) + { + bits = sign | ((unsigned int)(k + 150) << 23) | (m - (1U << 23)); + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + + r = -149 - k; + if (r >= 25) + { + bits = sign; + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + f = (unsigned int)vl_rne_u64((unsigned long long)m, r); + if (f == 0) + { + bits = sign; + } + else if (f == (1U << 23)) + { + bits = sign | (1U << 23); + } + else + { + bits = sign | f; + } + __builtin_memcpy(&x, &bits, sizeof x); + return x; +} + +/* + * Scale an x86 80-bit extended value by 2^n. The significand m is 64 bits + * wide with an explicit integer bit; se holds the sign (bit 15) and the + * exponent biased by 16383 (bits 14..0). A canonical normal has m in + * [2^63, 2^64) and p = ef - 16446 (value = m * 2^(ef - 16383 - 63)). A + * subnormal (ef == 0, m != 0) — and, defensively, an unnormal (ef > 0 with + * m < 2^63) — is interpreted as if the exponent were max(ef, 1) and + * normalizes exactly by shifting m left until its msb sits at bit 63. + * + * k >= 16321 overflows, k in [-16445, 16320] is an exact normal result + * (for the subnormal minimum, m = 2^63 with ef' = 1), and k <= -16446 is + * subnormal-or-zero with r = -16445 - k; r >= 65 always rounds to zero. + * The 64-bit significand makes r == 64 the deepest meaningful shift: only + * m above 2^63 then rounds up (to the smallest subnormal, 2^-16445). + */ +static inline long double +vl_scale2_ld(long double x, int n, int *overflowed) // NOLINT(bugprone-easily-swappable-parameters) +{ + struct + { + unsigned long long m; + unsigned short se; + } p; + unsigned long long m; + unsigned long long f; + unsigned short sign; + int ef; + int s; + int pwr; + int k; + int r; + + __builtin_memcpy(&p, &x, sizeof p); + m = p.m; + sign = (unsigned short)(p.se & 0x8000); + ef = p.se & 0x7fff; + + if (ef == 0x7fff) + { + return x; /* +-Inf and NaN pass through unchanged */ + } + if (m == 0) + { + return x; /* +-0 (and degenerate empty encodings) stay put */ + } + if (m < 0x8000000000000000ULL) + { + s = 63 - __builtin_clzll(m); + m <<= (63 - s); + pwr = (ef == 0 ? 1 : ef) + s - 16509; + } + else + { + pwr = ef - 16446; + } + + k = pwr + n; + if (k >= 16321) + { + if (overflowed != NULL) + { + *overflowed = 1; + } + p.se = (unsigned short)(sign | 0x7fff); + p.m = 0x8000000000000000ULL; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + if (k >= -16445) + { + p.se = (unsigned short)(sign | (unsigned short)(k + 16446)); + p.m = m; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + + r = -16445 - k; + if (r >= 65) + { + p.se = sign; + p.m = 0; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + f = vl_rne_u64(m, r); + p.se = sign; + p.m = f; + __builtin_memcpy(&x, &p, sizeof p); + return x; +} + +#endif /* VLIBC_MATH_IMPL_H */ diff --git a/src/math/modf.c b/src/math/modf.c new file mode 100644 index 0000000..8cad351 --- /dev/null +++ b/src/math/modf.c @@ -0,0 +1,215 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Split x into an integral part stored in *iptr and a fractional part + * returned, both carrying x's sign (C23 7.12.6.5): modf(-1.5, &i) puts + * -1.0 in i and returns -0.5. The integral part is x truncated toward + * zero. modf(±0, &i) stores ±0 and returns ±0; modf(±Inf, &i) stores + * ±Inf and returns ±0 with x's sign; modf(NaN, &i) stores the NaN and + * returns it; a subnormal |x| < 1 stores ±0 (x's sign) and returns x. + * + * The implementation is a bit-level truncation (the same "clear the + * fractional mantissa bits" shape as trunc.c): a value with |x| < 1 has + * integral part ±0; a value too large to carry a fraction (|x| >= 2^52, + * >= 2^23, >= 2^63 for the three formats) is its own integral part; the + * values in between have their low (fraction) bits masked off. The + * fractional part is then x - i, an exact subtraction — x and i share the + * significand's leading bits, so the difference is a representable + * multiple of the common unit in the last place. When the subtraction is + * an exact zero the fraction is ±0 with x's sign, which the hardware + * would otherwise always report as +0. + */ + +/* + * The double format: bit 63 the sign, bits 62..52 the exponent biased by + * 1023, bits 51..0 the fraction. A normal with exponent field e in + * [1023, 1074] has its low (1075 - e) fraction bits below the binary + * point; e >= 1075 means |x| >= 2^52, an integer already. + */ +double +modf(double x, double *iptr) +{ + unsigned long long bits; + unsigned long long sign; + unsigned long long ef; + unsigned long long frac; + unsigned long long ibits; + double i; + double f; + int shift; + + __builtin_memcpy(&bits, &x, sizeof bits); + sign = bits & (1ULL << 63); + ef = (bits >> 52) & 0x7ff; + frac = bits & 0xFFFFFFFFFFFFFULL; + + if (ef == 0x7ff) + { + *iptr = x; + if (frac != 0) + { + return x; /* NaN: both parts are the NaN */ + } + return sign != 0 ? -0.0 : 0.0; /* ±Inf: fraction ±0 of x's sign */ + } + if (ef == 0) + { + if (frac == 0) + { + *iptr = x; /* ±0 */ + return x; + } + *iptr = sign != 0 ? -0.0 : 0.0; /* subnormal: |x| < 1 */ + return x; + } + if (ef >= 1075) + { + *iptr = x; /* |x| >= 2^52: already integral */ + return sign != 0 ? -0.0 : 0.0; + } + if (ef < 1023) + { + *iptr = sign != 0 ? -0.0 : 0.0; /* 0 < |x| < 1 */ + return x; + } + + shift = 1075 - (int)ef; + ibits = bits & ~((1ULL << shift) - 1ULL); + __builtin_memcpy(&i, &ibits, sizeof i); + *iptr = i; + f = x - i; + if (f == 0.0) + { + return sign != 0 ? -0.0 : 0.0; + } + return f; +} + +/* + * The float format: bit 31 the sign, bits 30..23 the exponent biased by + * 127, bits 22..0 the fraction. For e in [127, 149] the low (150 - e) + * fraction bits are fractional; e >= 150 means |x| >= 2^23, integral. + */ +float +modff(float x, float *iptr) +{ + unsigned int bits; + unsigned int sign; + unsigned int ef; + unsigned int frac; + unsigned int ibits; + float i; + float f; + int shift; + + __builtin_memcpy(&bits, &x, sizeof bits); + sign = bits & (1U << 31); + ef = (bits >> 23) & 0xff; + frac = bits & 0x7FFFFFU; + + if (ef == 0xff) + { + *iptr = x; + if (frac != 0) + { + return x; + } + return sign != 0 ? -0.0f : 0.0f; + } + if (ef == 0) + { + if (frac == 0) + { + *iptr = x; + return x; + } + *iptr = sign != 0 ? -0.0f : 0.0f; + return x; + } + if (ef >= 150) + { + *iptr = x; + return sign != 0 ? -0.0f : 0.0f; + } + if (ef < 127) + { + *iptr = sign != 0 ? -0.0f : 0.0f; + return x; + } + + shift = 150 - (int)ef; + ibits = bits & ~((1U << shift) - 1U); + __builtin_memcpy(&i, &ibits, sizeof i); + *iptr = i; + f = x - i; + if (f == 0.0f) + { + return sign != 0 ? -0.0f : 0.0f; + } + return f; +} + +/* + * The x86 80-bit extended format: 64 significand bits m with an explicit + * integer bit and a sign/exponent word se biased by 16383. For e in + * [16383, 16445] the low (16446 - e) bits of m are fractional; e >= 16446 + * means |x| >= 2^63, integral. + */ +long double +modfl(long double x, long double *iptr) +{ + struct + { + unsigned long long m; + unsigned short se; + } p; + unsigned long long ibits; + long double i; + long double f; + int ef; + int shift; + + __builtin_memcpy(&p, &x, sizeof p); + ef = p.se & 0x7fff; + + if (ef == 0x7fff) + { + *iptr = x; + if (p.m != 0x8000000000000000ULL) + { + return x; /* NaN */ + } + return (p.se & 0x8000) != 0 ? -0.0L : 0.0L; /* ±Inf */ + } + if (p.m == 0) + { + *iptr = x; /* ±0 */ + return x; + } + if (ef >= 16446) + { + *iptr = x; + return (p.se & 0x8000) != 0 ? -0.0L : 0.0L; + } + if (ef < 16383) + { + *iptr = (p.se & 0x8000) != 0 ? -0.0L : 0.0L; /* subnormal: |x| < 1 */ + return x; + } + + shift = 16446 - ef; + ibits = p.m & ~((1ULL << shift) - 1ULL); + p.m = ibits; + __builtin_memcpy(&i, &p, sizeof i); + *iptr = i; + f = x - i; + if (f == 0.0L) + { + return (p.se & 0x8000) != 0 ? -0.0L : 0.0L; + } + return f; +} diff --git a/src/math/nearbyint.c b/src/math/nearbyint.c new file mode 100644 index 0000000..bff4630 --- /dev/null +++ b/src/math/nearbyint.c @@ -0,0 +1,248 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Integral value nearest to x in the current rounding direction (C23 + * 7.12.9.5), all three precisions, guaranteed never to raise the inexact + * exception. No exists in vlibc yet, so the only reachable + * rounding mode is the hardware default, round-to-nearest-even; the + * functions must still fold to correct results once lands, and + * the const attribute promises no exception and no errno path. + * + * GCC never folds the __builtin_nearbyint forms on this target (the + * roundsd expansion needs SSE4.1, absent from the default -march, so + * every __builtin_nearbyint call becomes an external nearbyint@PLT call + * at every optimization level), so each function below rounds directly + * on the IEEE 754 bit pattern with round-to-nearest-even: drop the low + * (frac-bits) of the significand; a dropped part above half an ulp steps + * the kept significand up, a dropped part of exactly half an ulp steps + * up only when the kept LSB is odd (ties to even). The step is a plain + * integer add that carries into the exponent field when the kept + * fraction is all ones, which renormalizes exactly; nothing here touches + * a rounding mode or raises an exception. nearbyint(±0) is ±0, + * nearbyint(-0.5) is -0.0 (tie to even zero), and ±Inf/NaN pass through. + */ + +/* + * The double format: bit 63 is the sign, bits 62..52 the exponent biased + * by 1023, bits 51..0 the fraction. A finite value with exponent e has + * fractional bits only when e < 1075; for e in [1023, 1074] exactly the + * low (1075 - e) bits of the fraction word are the fractional part. The + * significand is 53 bits wide (implicit 1 plus the 52-bit fraction), so + * in the e == 1023 binade the tie-even test looks at the implicit bit: + * the only half-way value there is 1.5, which rounds up to 2. + */ +static double +nearbyint_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 +nearbyint_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. A value m * 2^(e - 16446) has fractional bits + * only when e < 16446; 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 +nearbyint_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 nearbyint, for a float argument. + */ +float +nearbyintf(float x) +{ + return nearbyint_f(x); +} + +/* + * As nearbyint, for a double argument. + */ +double +nearbyint(double x) +{ + return nearbyint_d(x); +} + +/* + * As nearbyint, for a long double argument. + */ +long double +nearbyintl(long double x) +{ + return nearbyint_ld(x); +} diff --git a/src/math/remainder.c b/src/math/remainder.c new file mode 100644 index 0000000..711aee7 --- /dev/null +++ b/src/math/remainder.c @@ -0,0 +1,284 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#ifdef HAVE_CONFIG_H +#include +#endif + +/* + * The IEEE remainder x - n*y with n = x/y rounded to the nearest integer, + * ties to even (C23 7.12.10.2), all three precisions. + * + * The implementation mirrors the host glibc 2.44 remainder, which is a + * two-stage reduction: the exact fmod(x, y + y) first clears every whole + * multiple of 2y, leaving a residual below 2y (the doubling is finite by + * construction of the branch below), and the residual is then rounded to + * the nearest multiple of y -- at most one subtraction can still be wrong, + * so comparing the doubled residual against y decides the final n exactly, + * and each subtraction is exact (Sterbenz), making the result exact. The + * two stages run in x's magnitude and the sign of x is applied at the + * end, matching glibc. + * + * For |y| so large that y + y would overflow (|y| >= 2^1023 for double, + * >= 2^127 for float, >= 2^16383 for the 80-bit format) the fmod stage is + * skipped and the same rounding logic runs against y/2 -- with x below + * 2y at those magnitudes, at most two y subtractions can be needed, which + * the doubled-compare decision still performs exactly. Domain errors + * return a NaN (the classic indefinite pattern) and set errno to EDOM in + * the library build, matching the measured host behavior. + */ + +/* The double format: bit 63 the sign, bits 62..52 the exponent biased by + * 1023, bits 51..0 the fraction. */ +static double +remainder_d(double x, double y) +{ + const unsigned long long dbl_inf = 0x7FF0000000000000ULL; + unsigned long long xw; + unsigned long long yw; + unsigned long long hx; + unsigned long long hy; + unsigned long long sx; + double v; + double yh; + + __builtin_memcpy(&xw, &x, sizeof xw); + __builtin_memcpy(&yw, &y, sizeof yw); + sx = xw >> 63; + hx = xw & 0x7FFFFFFFFFFFFFFFULL; + hy = yw & 0x7FFFFFFFFFFFFFFFULL; + + v = fabs(x); + y = fabs(y); + if (hy < 0x7FE0000000000000ULL) + { + v = fmod(v, y + y); /* now v < 2y */ + if (v + v > y) + { + v -= y; + if (v + v >= y) + { + v -= y; + } + else if (v == 0.0) + { + v = 0.0; + } + } + } + else + { + /* x not finite or y a NaN */ + if (hx >= dbl_inf || hy > dbl_inf) + { +#ifdef HAVE_CONFIG_H + errno = EDOM; +#endif + xw = 0xFFF8000000000000ULL; + __builtin_memcpy(&x, &xw, sizeof xw); + return x; + } + yh = 0.5 * y; + if (v > yh) + { + v -= y; + if (v >= yh) + { + v -= y; + } + else if (v == 0.0) + { + v = 0.0; + } + } + } + + if (sx) + { + v = -v; + } + return v; +} + +/* The float format: bit 31 the sign, bits 30..23 the exponent biased by + * 127, bits 22..0 the fraction. */ +static float +remainder_f(float x, float y) +{ + const unsigned int flt_inf = 0x7F800000U; + unsigned int xw; + unsigned int yw; + unsigned int hx; + unsigned int hy; + unsigned int sx; + float v; + float yh; + + __builtin_memcpy(&xw, &x, sizeof xw); + __builtin_memcpy(&yw, &y, sizeof yw); + sx = xw >> 31; + hx = xw & 0x7FFFFFFFU; + hy = yw & 0x7FFFFFFFU; + + v = fabsf(x); + y = fabsf(y); + if (hy < 0x7F000000U) + { + v = fmodf(v, y + y); /* now v < 2y */ + if (v + v > y) + { + v -= y; + if (v + v >= y) + { + v -= y; + } + else if (v == 0.0f) + { + v = 0.0f; + } + } + } + else + { + /* x not finite or y a NaN */ + if (hx >= flt_inf || hy > flt_inf) + { +#ifdef HAVE_CONFIG_H + errno = EDOM; +#endif + xw = 0xFFC00000U; + __builtin_memcpy(&x, &xw, sizeof xw); + return x; + } + yh = 0.5f * y; + if (v > yh) + { + v -= y; + if (v >= yh) + { + v -= y; + } + else if (v == 0.0f) + { + v = 0.0f; + } + } + } + + if (sx) + { + v = -v; + } + return v; +} + +/* + * The x86 80-bit extended format: 64-bit significand m in bytes 0..7 (the + * integer bit is explicit) and a sign/exponent word se in bytes 8..9 with + * the sign in bit 15 and the exponent biased by 16383 in bits 14..0. + */ +struct remainder_ld_word +{ + unsigned long long m; + unsigned short se; +}; + +static long double +remainder_l(long double x, long double y) +{ + struct remainder_ld_word p; + struct remainder_ld_word q; + int sx; + int ex; + int ey; + long double v; + long double yh; + + __builtin_memcpy(&p, &x, sizeof p); + __builtin_memcpy(&q, &y, sizeof q); + sx = (int)(p.se >> 15); + ex = p.se & 0x7fff; + ey = q.se & 0x7fff; + + v = fabsl(x); + y = fabsl(y); + if (ey < 0x7ffe) + { + v = fmodl(v, y + y); /* now v < 2y */ + if (v + v > y) + { + v -= y; + if (v + v >= y) + { + v -= y; + } + else if (v == 0.0L) + { + v = 0.0L; + } + } + } + else + { + /* x not finite or y a NaN */ + if (ex == 0x7fff || (ey == 0x7fff && (q.m & 0x7FFFFFFFFFFFFFFFULL) != 0)) + { +#ifdef HAVE_CONFIG_H + errno = EDOM; +#endif + p.m = 0xC000000000000000ULL; + p.se = 0xFFFF; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + yh = 0.5L * y; + if (v > yh) + { + v -= y; + if (v >= yh) + { + v -= y; + } + else if (v == 0.0L) + { + v = 0.0L; + } + } + } + + if (sx) + { + v = -v; + } + return v; +} + +/* + * remainder for a float argument. + */ +float +remainderf(float x, float y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return remainder_f(x, y); +} + +/* + * remainder for a double argument. + */ +double +remainder(double x, double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return remainder_d(x, y); +} + +/* + * remainder for a long double argument. + */ +long double +remainderl(long double x, long double y) // NOLINT(bugprone-easily-swappable-parameters) +{ + return remainder_l(x, y); +} diff --git a/src/math/remquo.c b/src/math/remquo.c new file mode 100644 index 0000000..421f85d --- /dev/null +++ b/src/math/remquo.c @@ -0,0 +1,364 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * The IEEE remainder x - n*y with n = x/y rounded to nearest, ties to + * even (C23 7.12.10.2), plus the signed low bits of the integer quotient + * n stored through quo (C23 7.12.10.3), all three precisions. + * + * The structure mirrors the host glibc remainder/remquo implementation + * (glibc 2.44, measured byte-for-byte), because remquo's quotient is not + * the full n: glibc tracks only the three lowest quotient bits and still + * returns a remainder in the exact IEEE rounding, so the two functions + * share one reduce-then-round shape. x is first reduced by fmod(x, 8y) + * whenever |y| is small enough for 8y to be finite, which makes |x| < 8y + * and therefore bounds x/y below 8; the tracked quotient cquo collects + * the one 4y subtraction, the one 2y subtraction and the final rounding's + * up-to-two y subtractions, which together hold the three lowest bits of + * n. The final compare against y/2 (using doubled compares when y is so + * small that halving could lose precision) rounds x to the nearest + * multiple of y, ties to even via the natural 0.5 boundary, and each + * subtraction that fires is exact (Sterbenz), so the returned remainder + * is exact. + * + * The x86-64 host collapses the remainder and quotient signs the same + * way the code below does: the magnitude result is computed from |x| and + * |y| and then negated when x < 0, and *quo carries the sign of x/y. + * Domain errors -- |y| == 0, x not finite, y a NaN -- return a NaN (the + * classic indefinite pattern) and leave *quo untouched, as glibc does. + * Like glibc, remquo never touches errno. + */ + +/* The double format: bit 63 the sign, bits 62..52 the exponent biased by + * 1023, bits 51..0 the fraction. */ +static double +remquo_d(double x, double y, int *quo) +{ + const unsigned long long dbl_inf = 0x7FF0000000000000ULL; + unsigned long long xw; + unsigned long long yw; + unsigned long long hx; + unsigned long long hy; + unsigned long long sx; + unsigned long long sy; + double ax; + double ay; + double yh; + int cquo; + int qneg; + + __builtin_memcpy(&xw, &x, sizeof xw); + __builtin_memcpy(&yw, &y, sizeof yw); + sx = xw >> 63; + sy = yw >> 63; + qneg = (int)(sx ^ sy); + hx = xw & 0x7FFFFFFFFFFFFFFFULL; + hy = yw & 0x7FFFFFFFFFFFFFFFULL; + + if (hy == 0 || hx >= dbl_inf || hy > dbl_inf) + { + xw = 0xFFF8000000000000ULL; + __builtin_memcpy(&x, &xw, sizeof xw); + return x; + } + + if (hx == hy) + { + *quo = qneg ? -1 : 1; + return x * 0.0; + } + + ax = fabs(x); + ay = fabs(y); + if (hy <= 0x7FBFFFFFFFFFFFFFULL) + { + ax = fmod(ax, 8.0 * ay); /* now ax < 8 * ay */ + } + + cquo = 0; + if (hy <= 0x7FCFFFFFFFFFFFFFULL && ax >= 4.0 * ay) + { + ax -= 4.0 * ay; + cquo += 4; + } + if (hy <= 0x7FDFFFFFFFFFFFFFULL && ax >= 2.0 * ay) + { + ax -= 2.0 * ay; + cquo += 2; + } + + if (hy < 0x0020000000000000ULL) + { + if (ax + ax > ay) + { + ax -= ay; + cquo++; + if (ax + ax >= ay) + { + ax -= ay; + cquo++; + } + } + } + else + { + yh = 0.5 * ay; + if (ax > yh) + { + ax -= ay; + cquo++; + if (ax >= yh) + { + ax -= ay; + cquo++; + } + } + } + + *quo = qneg ? -cquo : cquo; + if (ax == 0.0) + { + ax = 0.0; + } + if (sx) + { + ax = -ax; + } + return ax; +} + +/* The float format: bit 31 the sign, bits 30..23 the exponent biased by + * 127, bits 22..0 the fraction. */ +static float +remquo_f(float x, float y, int *quo) +{ + const unsigned int flt_inf = 0x7F800000U; + unsigned int xw; + unsigned int yw; + unsigned int hx; + unsigned int hy; + unsigned int sx; + unsigned int sy; + float ax; + float ay; + float yh; + int cquo; + int qneg; + + __builtin_memcpy(&xw, &x, sizeof xw); + __builtin_memcpy(&yw, &y, sizeof yw); + sx = xw >> 31; + sy = yw >> 31; + qneg = (int)(sx ^ sy); + hx = xw & 0x7FFFFFFFU; + hy = yw & 0x7FFFFFFFU; + + if (hy == 0 || hx >= flt_inf || hy > flt_inf) + { + xw = 0xFFC00000U; + __builtin_memcpy(&x, &xw, sizeof xw); + return x; + } + + if (hx == hy) + { + *quo = qneg ? -1 : 1; + return x * 0.0f; + } + + ax = fabsf(x); + ay = fabsf(y); + if (hy <= 0x7DFFFFFFU) + { + ax = fmodf(ax, 8.0f * ay); /* now ax < 8 * ay */ + } + + cquo = 0; + if (hy <= 0x7E7FFFFFU && ax >= 4.0f * ay) + { + ax -= 4.0f * ay; + cquo += 4; + } + if (hy <= 0x7EFFFFFFU && ax >= 2.0f * ay) + { + ax -= 2.0f * ay; + cquo += 2; + } + + if (hy < 0x01000000U) + { + if (ax + ax > ay) + { + ax -= ay; + cquo++; + if (ax + ax >= ay) + { + ax -= ay; + cquo++; + } + } + } + else + { + yh = 0.5f * ay; + if (ax > yh) + { + ax -= ay; + cquo++; + if (ax >= yh) + { + ax -= ay; + cquo++; + } + } + } + + *quo = qneg ? -cquo : cquo; + if (ax == 0.0f) + { + ax = 0.0f; + } + if (sx) + { + ax = -ax; + } + return ax; +} + +/* + * The x86 80-bit extended format: 64-bit significand m in bytes 0..7 (the + * integer bit is explicit) and a sign/exponent word se in bytes 8..9 with + * the sign in bit 15 and the exponent biased by 16383 in bits 14..0. + */ +struct remquo_ld_word +{ + unsigned long long m; + unsigned short se; +}; + +static long double +remquo_l(long double x, long double y, int *quo) +{ + struct remquo_ld_word p; + struct remquo_ld_word q; + int ex; + int ey; + int qneg; + int sx; + int sy; + long double ax; + long double ay; + long double yh; + int cquo; + + __builtin_memcpy(&p, &x, sizeof p); + __builtin_memcpy(&q, &y, sizeof q); + sx = (int)(p.se >> 15); + sy = (int)(q.se >> 15); + qneg = sx ^ sy; + ex = p.se & 0x7fff; + ey = q.se & 0x7fff; + + if ((q.m == 0 && ey == 0) || ex == 0x7fff || + (ey == 0x7fff && (q.m & 0x7FFFFFFFFFFFFFFFULL) != 0)) + { + p.m = 0xC000000000000000ULL; + p.se = 0xFFFF; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + + if (ex == ey && p.m == q.m) + { + *quo = qneg ? -1 : 1; + return x * 0.0L; + } + + ax = fabsl(x); + ay = fabsl(y); + if (ey <= 0x7ffb) + { + ax = fmodl(ax, 8.0L * ay); /* now ax < 8 * ay */ + } + + cquo = 0; + if (ey <= 0x7ffc && ax >= 4.0L * ay) + { + ax -= 4.0L * ay; + cquo += 4; + } + if (ey <= 0x7ffd && ax >= 2.0L * ay) + { + ax -= 2.0L * ay; + cquo += 2; + } + + if (ey < 0x0002) + { + if (ax + ax > ay) + { + ax -= ay; + cquo++; + if (ax + ax >= ay) + { + ax -= ay; + cquo++; + } + } + } + else + { + yh = 0.5L * ay; + if (ax > yh) + { + ax -= ay; + cquo++; + if (ax >= yh) + { + ax -= ay; + cquo++; + } + } + } + + *quo = qneg ? -cquo : cquo; + if (ax == 0.0L) + { + ax = 0.0L; + } + if (sx) + { + ax = -ax; + } + return ax; +} + +/* + * remquo for a float argument. + */ +float +remquof(float x, float y, int *quo) // NOLINT(bugprone-easily-swappable-parameters) +{ + return remquo_f(x, y, quo); +} + +/* + * remquo for a double argument. + */ +double +remquo(double x, double y, int *quo) // NOLINT(bugprone-easily-swappable-parameters) +{ + return remquo_d(x, y, quo); +} + +/* + * remquo for a long double argument. + */ +long double +remquol(long double x, long double y, int *quo) // NOLINT(bugprone-easily-swappable-parameters) +{ + return remquo_l(x, y, quo); +} diff --git a/src/math/rint.c b/src/math/rint.c new file mode 100644 index 0000000..4abed0c --- /dev/null +++ b/src/math/rint.c @@ -0,0 +1,243 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * 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 + * 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); +} diff --git a/src/math/round.c b/src/math/round.c new file mode 100644 index 0000000..59cc863 --- /dev/null +++ b/src/math/round.c @@ -0,0 +1,226 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Integral value nearest to x, halfway cases rounded away from zero + * (C23 7.12.9.6), all three precisions. + * + * GCC never folds the __builtin_round forms on this target (no SSE4.1 in + * the default -march, so the roundsd expansion is unavailable and every + * __builtin_round call becomes an external round@PLT call at every + * optimization level), so each function below rounds directly on the + * IEEE 754 bit pattern. A negative argument with a nonzero dropped + * fraction must round away from zero, so the rounding is applied to the + * magnitude (sign bit untouched): drop the fractional bits, and when the + * dropped part is >= half an ulp step the kept significand up by one. + * The step is a plain integer add that carries into the exponent field + * when the kept fraction is all ones, which renormalizes exactly. + * round(±0) is ±0, round(-0.5) is -1.0 (half away from zero), and ±Inf + * and NaN pass through unchanged. Pure: no rounding mode consulted and no + * exception raised. + */ + +/* + * The double format: bit 63 is the sign, bits 62..52 the exponent biased + * by 1023, bits 51..0 the fraction. A finite value with exponent e has + * fractional bits only when e - 1023 < 52, i.e. e < 1075; for e in + * [1023, 1074] exactly the low (1075 - e) bits of the fraction word are + * the fractional part. + */ +static double +round_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); + + /* |x| >= 2^52 is already integral; Inf (e == 0x7ff) and NaN pass + * through unchanged as well. */ + if (e >= 1075) + { + return x; + } + + /* |x| < 1: round(±0) is ±0 and any other value with |x| >= 1/2 goes + * to ±1 while smaller magnitudes collapse to ±0, keeping the sign. */ + if (e < 1023) + { + if ((bits & ~(1ULL << 63)) == 0) + { + return x; + } + if (e == 1022) + { + 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) + { + 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. Fraction bits exist exactly when the + * exponent e is in [127, 149]; the low (150 - e) bits are fractional. + */ +static float +round_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) + { + 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) + { + 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) in bytes 0..7 and a sign/exponent word se in bytes 8..9, + * with the sign in bit 15 and the exponent (biased by 16383) in bits + * 14..0. A value m * 2^(e - 16446) has fractional bits only when + * e - 16383 < 63, i.e. e < 16446. The half-away step may overflow the + * 64-bit m when the kept significand is all ones; the carry then moves + * the value to the next binade (e + 1, integer bit alone). + */ +static long double +round_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) + { + 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) + { + 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 round, for a float argument. + */ +float +roundf(float x) +{ + return round_f(x); +} + +/* + * As round, for a double argument. + */ +double +round(double x) +{ + return round_d(x); +} + +/* + * As round, for a long double argument. + */ +long double +roundl(long double x) +{ + return round_ld(x); +} diff --git a/src/math/scalbn.c b/src/math/scalbn.c new file mode 100644 index 0000000..bee1486 --- /dev/null +++ b/src/math/scalbn.c @@ -0,0 +1,121 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#include "math_impl.h" + +/* + * x * FLT_RADIX^n with FLT_RADIX 2 (C23 7.12.6.6), all three precisions + * and both exponent-argument types. scalbn takes an int n, scalbln a long + * n; otherwise the two families are semantically identical to ldexp (they + * share the exact same scale cores and overflow-to-+-Inf-with-ERANGE + * behavior, and subnormal inputs are handled exactly like subnormal + * outputs). The scalbln functions clamp the long exponent to +-20000 + * first: any magnitude beyond that saturates every result to +-Inf or +-0 + * in all three formats, and clamping keeps the arithmetic inside int range + * with no shift by a huge count. + */ + +/* + * As ldexp (see ldexp.c) for a double x and an int n. + */ +double +scalbn(double x, int n) +{ + int overflowed = 0; + double r = vl_scale2_d(x, n, &overflowed); + +#ifdef HAVE_CONFIG_H + if (overflowed) + { + errno = ERANGE; + } +#endif + return r; +} + +float +scalbnf(float x, int n) +{ + int overflowed = 0; + float r = vl_scale2_f(x, n, &overflowed); + +#ifdef HAVE_CONFIG_H + if (overflowed) + { + errno = ERANGE; + } +#endif + return r; +} + +long double +scalbnl(long double x, int n) +{ + int overflowed = 0; + long double r = vl_scale2_ld(x, n, &overflowed); + +#ifdef HAVE_CONFIG_H + if (overflowed) + { + errno = ERANGE; + } +#endif + return r; +} + +/* + * The scalbln family: as scalbn with the exponent given as a long. + * n is first clamped into [-20000, 20000]; anything beyond saturates every + * precision's range (the largest meaningful long-double exponent is below + * 16446 in magnitude), so no precision is lost by the clamp. + */ +double +scalbln(double x, long n) +{ + if (n > 20000) + { + n = 20000; + } + else if (n < -20000) + { + n = -20000; + } + return scalbn(x, (int)n); +} + +float +scalblnf(float x, long n) +{ + if (n > 20000) + { + n = 20000; + } + else if (n < -20000) + { + n = -20000; + } + return scalbnf(x, (int)n); +} + +long double +scalblnl(long double x, long n) +{ + if (n > 20000) + { + n = 20000; + } + else if (n < -20000) + { + n = -20000; + } + return scalbnl(x, (int)n); +} diff --git a/src/math/trunc.c b/src/math/trunc.c new file mode 100644 index 0000000..acc3e93 --- /dev/null +++ b/src/math/trunc.c @@ -0,0 +1,127 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +/* + * Integral value nearest to x in the direction of zero (C23 7.12.9.3), + * all three precisions. + * + * GCC does not fold the __builtin_trunc forms on this target (see + * floor.c), so each function below simply clears the fractional mantissa + * bits on the IEEE 754 bit pattern. Rounding toward zero never rounds, so + * no rounding mode or exception is involved; trunc(-0.3) is -0.0 (the + * sign bit survives), and ±Inf/NaN pass through unchanged. + */ + +static double +trunc_d(double x) +{ + unsigned long long bits; + int e; + + __builtin_memcpy(&bits, &x, sizeof bits); + e = (int)((bits >> 52) & 0x7ff); + + /* |x| >= 2^52 is already integral; Inf (e == 0x7ff) and NaN pass + * through unchanged as well. */ + if (e >= 1075) + { + return x; + } + if (e < 1023) + { + /* |x| < 1: truncation is ±0, keeping the sign bit. */ + bits &= 1ULL << 63; + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + + bits &= ~((1ULL << (1075 - e)) - 1ULL); + __builtin_memcpy(&x, &bits, sizeof x); + return x; +} + +static float +trunc_f(float x) +{ + unsigned int bits; + int e; + + __builtin_memcpy(&bits, &x, sizeof bits); + e = (int)((bits >> 23) & 0xff); + + if (e >= 150) + { + return x; + } + if (e < 127) + { + bits &= 1U << 31; + __builtin_memcpy(&x, &bits, sizeof x); + return x; + } + + bits &= ~((1U << (150 - e)) - 1U); + __builtin_memcpy(&x, &bits, sizeof x); + return x; +} + +static long double +trunc_ld(long double x) +{ + struct + { + unsigned long long m; + unsigned short se; + } p; + int e; + + __builtin_memcpy(&p, &x, sizeof p); + e = p.se & 0x7fff; + + if (e >= 16446) + { + return x; + } + if (e < 16383) + { + /* |x| < 1: truncation is ±0, keeping the sign bit. */ + p.m = 0; + p.se &= 0x8000; + __builtin_memcpy(&x, &p, sizeof p); + return x; + } + + p.m &= ~((1ULL << (63 - (e - 16383))) - 1ULL); + __builtin_memcpy(&x, &p, sizeof p); + return x; +} + +/* + * As trunc, for a float argument. + */ +float +truncf(float x) +{ + return trunc_f(x); +} + +/* + * As trunc, for a double argument. + */ +double +trunc(double x) +{ + return trunc_d(x); +} + +/* + * As trunc, for a long double argument. + */ +long double +truncl(long double x) +{ + return trunc_ld(x); +} diff --git a/tests/test_math_basic.c b/tests/test_math_basic.c new file mode 100644 index 0000000..cc3c1b4 --- /dev/null +++ b/tests/test_math_basic.c @@ -0,0 +1,342 @@ +/* + * vlibc — math.h todo-39 consolidated final test: slices S1..S4. + * + * One pass over every implemented todo-39 family, in all three + * precisions, as a final smoke of the whole arithmetic slice: + * + * S1 fabs/fabsf/fabsl, copysign, floor, ceil, trunc + * S2 round, rint, nearbyint, lrint, llrint, ilogb, logb + * S3 frexp, ldexp, modf, scalbn, scalbln + * S4 fmin, fmax, fdim, lround, llround, fmod, remainder, remquo + * + * Expected values follow IEEE 754-2008 / C23 semantics by hand and, for + * the S4 families, the byte-exact host glibc 2.44 behavior that the + * differential probe checks (the remquo quotient rows, in particular, are + * glibc's tracked-low-bits values, e.g. remquo(100.0, 3.0) reports 1, the + * low quotient bits, not 33). Inputs travel through volatile objects so + * the compiler cannot constant-fold the calls away; diagnostics go + * through raw SYS_write and no host headers are included. + */ + +#include + +#include "../src/internal/syscall.h" + +/* Write a NUL-terminated string to fd via the raw syscall layer. */ +static __attribute__((optimize("no-tree-loop-distribute-patterns"))) void +say(int fd, const char *s) +{ + long n = 0; + + while (s[n] != '\0') + { + n++; + } + __syscall3(SYS_write, fd, (long)s, n); +} + +static int failures; + +static void +check(int ok, const char *msg) +{ + if (ok) + { + say(1, "ok "); + } + else + { + say(1, "FAIL "); + failures++; + } + say(1, msg); + say(1, "\n"); +} + +/* S1: fabs, copysign, floor, ceil, trunc. */ +static void +s1_scenario(void) +{ + volatile double negzero = -0.0; + volatile double nan = NAN; + volatile double big = 1e300; + volatile float nanf = NAN; + volatile long double nanl = NAN; + + check(fabs(-3.5) == 3.5, "fabs(-3.5) == 3.5"); + check(fabs(negzero) == 0.0 && !signbit(fabs(negzero)), "fabs(-0.0) == +0.0"); + check(fabs(nan) != nan, "fabs(NaN) is NaN"); + check(copysign(3.0, -1.0) == -3.0, "copysign(3,-1) == -3"); + check(copysign(-3.0, 1.0) == 3.0, "copysign(-3,1) == 3"); + check(copysign(negzero, 1.0) == 0.0 && !signbit(copysign(negzero, 1.0)), + "copysign(-0,+1) == +0"); + check(floor(3.7) == 3.0, "floor(3.7) == 3"); + check(floor(-3.7) == -4.0, "floor(-3.7) == -4"); + check(floor(negzero) == 0.0 && signbit(floor(negzero)), "floor(-0.0) == -0.0"); + check(ceil(3.2) == 4.0, "ceil(3.2) == 4"); + check(ceil(-3.2) == -3.0, "ceil(-3.2) == -3"); + check(trunc(3.7) == 3.0, "trunc(3.7) == 3"); + check(trunc(-3.7) == -3.0, "trunc(-3.7) == -3"); + check(fabsf(-2.5f) == 2.5f, "fabsf(-2.5f) == 2.5f"); + check(fabsf(nanf) != nanf, "fabsf(NaN) is NaN"); + check(copysignf(1.0f, -0.0f) == -1.0f, "copysignf(1,-0) == -1"); + check(floorf(2.9f) == 2.0f, "floorf(2.9f) == 2"); + check(ceilf(-2.1f) == -2.0f, "ceilf(-2.1f) == -2"); + check(truncf(2.9f) == 2.0f, "truncf(2.9f) == 2"); + check(fabsl(-3.5L) == 3.5L, "fabsl(-3.5L) == 3.5L"); + check(fabsl(nanl) != nanl, "fabsl(NaN) is NaN"); + check(copysignl(-4.0L, -2.0L) == -4.0L, "copysignl(-4,-2) == -4"); + check(floorl(-3.7L) == -4.0L, "floorl(-3.7L) == -4"); + check(ceill(3.2L) == 4.0L, "ceill(3.2L) == 4"); + check(truncl(-3.7L) == -3.0L, "truncl(-3.7L) == -3"); + check(floor(big) == big, "floor(1e300) == 1e300"); +} + +/* S2: round, rint, nearbyint, lrint, llrint, ilogb, logb. */ +static void +s2_scenario(void) +{ + volatile double negzero = -0.0; + volatile double nan = NAN; + volatile double sub = 0x1p-1074; + volatile double two40 = 0x1p40; + + check(round(0.5) == 1.0, "round(0.5) == 1"); + check(round(-0.5) == -1.0, "round(-0.5) == -1"); + check(round(2.5) == 3.0, "round(2.5) == 3"); + check(round(negzero) == 0.0 && signbit(round(negzero)), "round(-0.0) == -0.0"); + check(round(sub) == 0.0 && !signbit(round(sub)), "round(5e-324) == +0"); + check(isnan(round(nan)), "round(NaN) is NaN"); + check(rint(0.5) == 0.0, "rint(0.5) == 0 (ties to even)"); + check(rint(2.5) == 2.0, "rint(2.5) == 2 (ties to even)"); + check(rint(negzero) == 0.0 && signbit(rint(negzero)), "rint(-0.0) == -0.0"); + check(nearbyint(2.5) == 2.0, "nearbyint(2.5) == 2 (ties to even)"); + check(nearbyint(1.5) == 2.0, "nearbyint(1.5) == 2"); + check(lrint(2.5) == 2L, "lrint(2.5) == 2"); + check(lrint(two40) == 1099511627776L, "lrint(2^40) == 2^40"); + check(llrint(-2.5) == -2LL, "llrint(-2.5) == -2"); + check(llrint(1.5) == 2LL, "llrint(1.5) == 2"); + check(ilogb(8.0) == 3, "ilogb(8) == 3"); + check(ilogb(0.5) == -1, "ilogb(0.5) == -1"); + check(ilogb(sub) == -1074, "ilogb(min subnormal) == -1074"); + check(ilogb(0.0) == FP_ILOGB0, "ilogb(0.0) == FP_ILOGB0"); + check(logb(8.0) == 3.0, "logb(8) == 3"); + check(logb(0.5) == -1.0, "logb(0.5) == -1"); + check(logb(0.0) == -HUGE_VAL, "logb(0.0) == -Inf"); + check(roundf(-0.5f) == -1.0f, "roundf(-0.5f) == -1"); + check(rintf(2.5f) == 2.0f, "rintf(2.5f) == 2 (ties to even)"); + check(nearbyintf(-2.5f) == -2.0f, "nearbyintf(-2.5f) == -2"); + check(lrintf(1.5f) == 2L, "lrintf(1.5f) == 2"); + check(llrintf(2.5f) == 2LL, "llrintf(2.5f) == 2"); + check(ilogbf(0x1p-149f) == -149, "ilogbf(min subnormal) == -149"); + check(logbf(8.0f) == 3.0f, "logbf(8) == 3"); + check(roundl(2.5L) == 3.0L, "roundl(2.5L) == 3"); + check(rintl(-0.5L) == 0.0L && signbit(rintl(-0.5L)), "rintl(-0.5L) == -0.0L"); + check(nearbyintl(1.5L) == 2.0L, "nearbyintl(1.5L) == 2"); + check(lrintl(2.5L) == 2L, "lrintl(2.5L) == 2"); + check(llrintl(-2.5L) == -2LL, "llrintl(-2.5L) == -2"); + check(ilogbl(0x1p-16445L) == -16445, "ilogbl(min subnormal) == -16445"); + check(logbl(0.5L) == -1.0L, "logbl(0.5L) == -1"); +} + +/* S3: frexp, ldexp, modf, scalbn, scalbln. */ +static void +s3_scenario(void) +{ + volatile double negzero = -0.0; + int e; + double ip; + double fp; + float ipf; + float fpf; + long double ipl; + long double fpl; + + check(frexp(8.0, &e) == 0.5 && e == 4, "frexp(8) == (0.5, 4)"); + check(frexp(-6.0, &e) == -0.75 && e == 3, "frexp(-6) == (-0.75, 3)"); + check(frexp(0.0, &e) == 0.0 && e == 0, "frexp(0) == (0, 0)"); + check(frexp(negzero, &e) == 0.0 && e == 0, "frexp(-0) == (-0, 0)"); + check(ldexp(1.0, 10) == 1024.0, "ldexp(1,10) == 1024"); + check(ldexp(1.0, -1075) == 0.0, "ldexp(1,-1075) == +0"); + check(scalbn(1.0, 10) == 1024.0, "scalbn(1,10) == 1024"); + check(scalbln(1.0, 10L) == 1024.0, "scalbln(1,10) == 1024"); + fp = modf(3.75, &ip); + check(ip == 3.0 && fp == 0.75, "modf(3.75) == (3, 0.75)"); + fp = modf(-3.75, &ip); + check(ip == -3.0 && fp == -0.75, "modf(-3.75) == (-3, -0.75)"); + check(frexpf(8.0f, &e) == 0.5f && e == 4, "frexpf(8) == (0.5, 4)"); + check(ldexpf(1.0f, 10) == 1024.0f, "ldexpf(1,10) == 1024"); + check(scalbnf(1.0f, -10) == 0.0009765625f, "scalbnf(1,-10) == 2^-10"); + fpf = modff(3.75f, &ipf); + check(ipf == 3.0f && fpf == 0.75f, "modff(3.75f) == (3, 0.75)"); + check(frexpl(8.0L, &e) == 0.5L && e == 4, "frexpl(8) == (0.5, 4)"); + check(ldexpl(1.0L, 10) == 1024.0L, "ldexpl(1,10) == 1024"); + check(scalbnl(1.0L, 10) == 1024.0L, "scalbnl(1,10) == 1024"); + fpl = modfl(3.75L, &ipl); + check(ipl == 3.0L && fpl == 0.75L, "modfl(3.75L) == (3, 0.75)"); +} + +/* S4a: fmin, fmax, fdim, lround, llround. */ +static void +s4_minmax_scenario(void) +{ + volatile double negzero = -0.0; + volatile double nan = NAN; + volatile double inf = HUGE_VAL; + volatile float nanf = NAN; + volatile long double nanl = NAN; + volatile double big = 1e300; + volatile double huge = 1.7e308; + volatile float bigf = 1e30f; + + check(fmin(3.0, 5.0) == 3.0, "fmin(3,5) == 3"); + check(fmin(5.0, 3.0) == 3.0, "fmin(5,3) == 3"); + check(fmin(3.0, 5.0) == 3.0, "fmin(3,5) == 3"); + check(fmin(negzero, 0.0) == 0.0 && signbit(fmin(negzero, 0.0)), "fmin(-0,+0) == -0"); + check(fmin(0.0, negzero) == 0.0 && signbit(fmin(0.0, negzero)), "fmin(+0,-0) == -0"); + check(fmin(nan, 5.0) == 5.0, "fmin(NaN,5) == 5"); + check(fmin(5.0, nan) == 5.0, "fmin(5,NaN) == 5"); + check(isnan(fmin(nan, nan)), "fmin(NaN,NaN) is NaN"); + check(fmax(3.0, 5.0) == 5.0, "fmax(3,5) == 5"); + check(fmax(5.0, 3.0) == 5.0, "fmax(5,3) == 5"); + check(fmax(negzero, 0.0) == 0.0 && !signbit(fmax(negzero, 0.0)), "fmax(-0,+0) == +0"); + check(fmax(0.0, negzero) == 0.0 && !signbit(fmax(0.0, negzero)), "fmax(+0,-0) == +0"); + check(fmax(nan, 5.0) == 5.0, "fmax(NaN,5) == 5"); + check(fmax(5.0, nan) == 5.0, "fmax(5,NaN) == 5"); + check(fmin(-inf, -5.0) == -inf, "fmin(-Inf,-5) == -Inf"); + check(fmax(inf, 5.0) == inf, "fmax(+Inf,5) == +Inf"); + check(fminf(3.0f, 5.0f) == 3.0f, "fminf(3,5) == 3"); + check(fminf(nanf, 5.0f) == 5.0f, "fminf(NaN,5) == 5"); + check(fmaxf(3.0f, 5.0f) == 5.0f, "fmaxf(3,5) == 5"); + check(fminl(3.0L, 5.0L) == 3.0L, "fminl(3,5) == 3"); + check(fminl(nanl, 5.0L) == 5.0L, "fminl(NaN,5) == 5"); + check(fmaxl(3.0L, 5.0L) == 5.0L, "fmaxl(3,5) == 5"); + check(fmaxl(nanl, 5.0L) == 5.0L, "fmaxl(NaN,5) == 5"); + + check(fdim(5.0, 3.0) == 2.0, "fdim(5,3) == 2"); + check(fdim(3.0, 5.0) == 0.0 && !signbit(fdim(3.0, 5.0)), "fdim(3,5) == +0"); + check(fdim(3.0, 3.0) == 0.0, "fdim(3,3) == +0"); + check(fdim(inf, -inf) == inf, "fdim(+Inf,-Inf) == +Inf"); + check(isnan(fdim(nan, 5.0)), "fdim(NaN,5) is NaN"); + check(isnan(fdim(5.0, nan)), "fdim(5,NaN) is NaN"); + check(fdim(big, -big) == 2.0 * big, "fdim(1e300,-1e300) == 2e300"); + check(fdim(huge, -huge) == inf, "fdim(1.7e308,-1.7e308) == +Inf"); + check(fdimf(5.0f, 3.0f) == 2.0f, "fdimf(5,3) == 2"); + check(fdimf(3.0f, 5.0f) == 0.0f, "fdimf(3,5) == +0"); + check(fdiml(5.0L, 3.0L) == 2.0L, "fdiml(5,3) == 2"); + check(fdiml(3.0L, 5.0L) == 0.0L, "fdiml(3,5) == +0"); + + check(lround(0.5) == 1L, "lround(0.5) == 1"); + check(lround(-0.5) == -1L, "lround(-0.5) == -1"); + check(lround(2.5) == 3L, "lround(2.5) == 3"); + check(lround(-2.5) == -3L, "lround(-2.5) == -3"); + check(lround(big) == LONG_MIN, "lround(1e300) == LONG_MIN"); + check(lround(-big) == LONG_MIN, "lround(-1e300) == LONG_MIN"); + check(lroundf(2.5f) == 3L, "lroundf(2.5f) == 3"); + check(lroundf(bigf) == LONG_MIN, "lroundf(1e30f) == LONG_MIN"); + check(lroundl(2.5L) == 3L, "lroundl(2.5L) == 3"); + check(llround(0.5) == 1LL, "llround(0.5) == 1"); + check(llround(-0.5) == -1LL, "llround(-0.5) == -1"); + check(llround(2.5) == 3LL, "llround(2.5) == 3"); + check(llround(big) == LLONG_MIN, "llround(1e300) == LLONG_MIN"); + check(llround(nan) == LLONG_MIN, "llround(NaN) == LLONG_MIN"); + check(llroundf(2.5f) == 3LL, "llroundf(2.5f) == 3"); + check(llroundl(2.5L) == 3LL, "llroundl(2.5L) == 3"); +} + +/* S4b: fmod, remainder, remquo (glibc-verified rows). */ +static void +s4_rem_scenario(void) +{ + volatile double negzero = -0.0; + volatile double nan = NAN; + volatile double inf = HUGE_VAL; + volatile float nanf = NAN; + volatile long double nanl = NAN; + volatile double a100 = 100.0; + volatile double a34 = 34.0; + volatile double a40 = 40.0; + volatile double a3 = 3.0; + volatile double a5 = 5.0; + volatile float f100 = 100.0f; + volatile long double l100 = 100.0L; + volatile long double l40 = 40.0L; + int q; + + check(fmod(5.0, 3.0) == 2.0, "fmod(5,3) == 2"); + check(fmod(-5.0, 3.0) == -2.0, "fmod(-5,3) == -2"); + check(fmod(5.0, -3.0) == 2.0, "fmod(5,-3) == 2"); + check(fmod(-5.0, -3.0) == -2.0, "fmod(-5,-3) == -2"); + check(fmod(6.0, 3.0) == 0.0 && !signbit(fmod(6.0, 3.0)), "fmod(6,3) == +0"); + check(fmod(-6.0, 3.0) == 0.0 && signbit(fmod(-6.0, 3.0)), "fmod(-6,3) == -0"); + check(fmod(3.0, 5.0) == 3.0, "fmod(3,5) == 3"); + check(fmod(negzero, 5.0) == 0.0 && signbit(fmod(negzero, 5.0)), "fmod(-0,5) == -0"); + check(fmod(5.0, inf) == 5.0, "fmod(5,+Inf) == 5"); + check(isnan(fmod(5.0, 0.0)), "fmod(5,0) is NaN"); + check(isnan(fmod(inf, 3.0)), "fmod(+Inf,3) is NaN"); + check(isnan(fmod(nan, 3.0)), "fmod(NaN,3) is NaN"); + check(fmodf(5.0f, 3.0f) == 2.0f, "fmodf(5,3) == 2"); + check(fmodf(-5.0f, 3.0f) == -2.0f, "fmodf(-5,3) == -2"); + check(fmodf(5.0f, 0.0f) != 0 && !(fmodf(5.0f, 0.0f) == fmodf(5.0f, 0.0f)), "fmodf(5,0) is NaN"); + check(isnan(fmodf(nanf, 3.0f)), "fmodf(NaN,3) is NaN"); + check(fmodl(5.0L, 3.0L) == 2.0L, "fmodl(5,3) == 2"); + check(fmodl(-5.0L, 3.0L) == -2.0L, "fmodl(-5,3) == -2"); + check(fmodl(5.0L, 0.0L) != fmodl(5.0L, 0.0L), "fmodl(5,0) is NaN"); + check(isnan(fmodl(nanl, 3.0L)), "fmodl(NaN,3) is NaN"); + + check(remainder(5.0, 3.0) == -1.0, "remainder(5,3) == -1"); + check(remainder(5.0, 2.0) == 1.0, "remainder(5,2) == 1"); + check(remainder(6.0, 4.0) == -2.0, "remainder(6,4) == -2"); + check(remainder(7.0, 4.0) == -1.0, "remainder(7,4) == -1"); + check(remainder(5.0, 5.0) == 0.0, "remainder(5,5) == +0"); + check(remainder(-5.0, 3.0) == 1.0, "remainder(-5,3) == 1"); + check(remainder(5.0, -3.0) == -1.0, "remainder(5,-3) == -1"); + check(isnan(remainder(5.0, 0.0)), "remainder(5,0) is NaN"); + check(isnan(remainder(inf, 3.0)), "remainder(+Inf,3) is NaN"); + check(isnan(remainder(nan, 3.0)), "remainder(NaN,3) is NaN"); + check(remainderf(5.0f, 3.0f) == -1.0f, "remainderf(5,3) == -1"); + check(isnan(remainderf(5.0f, 0.0f)), "remainderf(5,0) is NaN"); + check(remainderl(5.0L, 3.0L) == -1.0L, "remainderl(5,3) == -1"); + check(isnan(remainderl(5.0L, 0.0L)), "remainderl(5,0) is NaN"); + check(remainderl(-5.0L, 3.0L) == 1.0L, "remainderl(-5,3) == 1"); + + check(remquo(5.0, 3.0, &q) == -1.0 && q == 2, "remquo(5,3): r=-1 q=2"); + check(remquo(-5.0, 3.0, &q) == 1.0 && q == -2, "remquo(-5,3): r=1 q=-2"); + check(remquo(5.0, -3.0, &q) == -1.0 && q == -2, "remquo(5,-3): r=-1 q=-2"); + check(remquo(a100, a3, &q) == 1.0 && q == 1, "remquo(100,3): r=1 q=1"); + check(remquo(a34, a3, &q) == 1.0 && q == 3, "remquo(34,3): r=1 q=3"); + check(remquo(a40, a3, &q) == 1.0 && q == 5, "remquo(40,3): r=1 q=5"); + check(remquo(7.0, 2.0, &q) == -1.0 && q == 4, "remquo(7,2): r=-1 q=4"); + check(remquo(6.0, 3.0, &q) == 0.0 && q == 2, "remquo(6,3): r=0 q=2"); + check(remquo(2.0, 3.0, &q) == -1.0 && q == 1, "remquo(2,3): r=-1 q=1"); + check(remquo(1.0, 3.0, &q) == 1.0 && q == 0, "remquo(1,3): r=1 q=0"); + check(remquo(a5, a5, &q) == 0.0 && q == 1, "remquo(5,5): r=0 q=1"); + check(isnan(remquo(a5, 0.0, &q)), "remquo(5,0) is NaN"); + check(isnan(remquo(inf, a3, &q)), "remquo(+Inf,3) is NaN"); + check(remquof(f100, 3.0f, &q) == 1.0f && q == 1, "remquof(100,3): q=1"); + check(remquof(5.0f, 3.0f, &q) == -1.0f && q == 2, "remquof(5,3): q=2"); + check(remquol(l100, 3.0L, &q) == 1.0L && q == 1, "remquol(100,3): q=1"); + check(remquol(5.0L, 3.0L, &q) == -1.0L && q == 2, "remquol(5,3): q=2"); + check(remquol(l40, 3.0L, &q) == 1.0L && q == 5, "remquol(40,3): q=5"); + check(isnan(remquol(nanl, 3.0L, &q)), "remquol(NaN,3) is NaN"); +} + +int +main(void) +{ + s1_scenario(); + s2_scenario(); + s3_scenario(); + s4_minmax_scenario(); + s4_rem_scenario(); + + if (failures == 0) + { + say(1, "all math todo-39 (S1-S4) tests passed\n"); + } + else + { + say(1, "FAILURES\n"); + } + return failures == 0 ? 0 : 1; +} diff --git a/tests/test_math_s1.c b/tests/test_math_s1.c new file mode 100644 index 0000000..ed3a0d1 --- /dev/null +++ b/tests/test_math_s1.c @@ -0,0 +1,264 @@ +/* + * vlibc — math.h slice S1 test: fabs/copysign/floor/ceil/trunc (todo 39). + * + * Exercises every implemented family in all three precisions plus the + * classification macros and the generic dispatch: + * + * 1. fabs/fabsf/fabsl of ±1.5 and of ±0.0 (signed-zero flattening); + * 2. copysign/copysignf/copysignl sign transfer, including from ±0.0 + * and onto a NaN; + * 3. floor across the sign boundary: floor(-0.5) == -1.0, floor(+0.5) + * == +0.0, integers and ±0 pass through, and the subnormal + * floor(5e-324) == 0.0 / floor(-5e-324) == -1.0; + * 4. ceil with the real negative zero: ceil(-0.5) == -0.0 with signbit + * set, ceil(+0.5) == 1.0; + * 5. trunc toward zero: trunc(-1.7) == -1.0, trunc(1.7) == 1.0, + * trunc(-0.7) == -0.0 (signbit set); + * 6. NaN and ±Inf passthrough for all three rounding functions and the + * boundary mantissa-clearing paths (values whose integer part fills + * the whole significand minus one bit); + * 7. classification macros (fpclassify/isnan/isinf/isfinite/isnormal/ + * signbit) over the FP_* constants. + * + * Values that must survive exactly are fed through volatile objects so the + * compiler cannot constant-fold the call away and the real library + * functions run. No host headers are included (-Iinclude shadows GCC's); + * diagnostics go through raw SYS_write, and errno is never read (vlibc's + * errno slot collides with the host TCB in this host-linked standalone + * binary). + */ + +#include +#include + +#include "../src/internal/syscall.h" + +/* Write a NUL-terminated string to fd via the raw syscall layer. The + * optimize attribute keeps GCC from lowering the length loop into a + * strlen call, which would leave a vlibc-owned symbol undefined in this + * host-linked standalone binary (house idiom, see src/string). */ +static __attribute__((optimize("no-tree-loop-distribute-patterns"))) void +say(int fd, const char *s) +{ + long n = 0; + + while (s[n] != '\0') + { + n++; + } + __syscall3(SYS_write, fd, (long)s, n); +} + +static int failures; + +static void +check(int ok, const char *msg) +{ + if (ok) + { + say(1, "ok "); + } + else + { + say(1, "FAIL "); + failures++; + } + say(1, msg); + say(1, "\n"); +} + +/* 1. fabs family: magnitude of positive and negative inputs. */ +static void +fabs_scenario(void) +{ + volatile float xf = -1.5f; + volatile double xd = -1.5; + volatile long double xl = -1.5L; + volatile float zf = -0.0f; + volatile double zd = -0.0; + + check(fabsf(xf) == 1.5f, "fabsf(-1.5f) == 1.5f"); + check(fabs(xd) == 1.5, "fabs(-1.5) == 1.5"); + check(fabsl(xl) == 1.5L, "fabsl(-1.5L) == 1.5L"); + check(fabsf(1.5f) == 1.5f, "fabsf(+1.5f) == 1.5f"); + check(fabs(1.5) == 1.5, "fabs(+1.5) == 1.5"); + check(fabsl(1.5L) == 1.5L, "fabsl(+1.5L) == 1.5L"); + check(fabsf(zf) == 0.0f && !signbit(fabsf(zf)), "fabsf(-0.0f) is +0.0"); + check(fabs(zd) == 0.0 && !signbit(fabs(zd)), "fabs(-0.0) is +0.0"); +} + +/* 2. copysign family: sign follows y, including y == ±0.0 and NaN x. */ +static void +copysign_scenario(void) +{ + volatile double negzero = -0.0; + volatile double poszero = 0.0; + volatile double nan = NAN; + + check(copysign(1.0, negzero) == -1.0, "copysign(1.0, -0.0) == -1.0"); + check(copysign(1.0, poszero) == 1.0, "copysign(1.0, +0.0) == +1.0"); + check(copysign(-1.0, poszero) == 1.0, "copysign(-1.0, +0.0) == +1.0"); + check(signbit(copysign(1.0, negzero)), "copysign(1.0, -0.0) has signbit set"); + check(copysignf(-1.0f, -0.0f) == -1.0f, "copysignf(-1.0f, -0.0f) == -1.0f"); + check(copysignl(-1.0L, 0.0L) == 1.0L, "copysignl(-1.0L, +0.0L) == +1.0L"); + check(signbit(copysign(nan, negzero)), "copysign(NaN, -0.0) carries the sign"); + check(!signbit(copysign(nan, poszero)), "copysign(NaN, +0.0) clears the sign"); +} + +/* 3. floor: toward -Inf, across the sign boundary and the subnormal + * range. */ +static void +floor_scenario(void) +{ + volatile double half = 0.5; + volatile long double neg_half = -0.5L; + volatile double two = 2.0; + volatile float twof = 2.0f; + volatile double negzero = -0.0; + volatile double sub = 5e-324; /* the smallest subnormal */ + volatile double big = 2251799813685248.0; /* 2^51 */ + volatile double big_half = 2251799813685248.5; /* 2^51 + 0.5 */ + + check(floor(-0.5) == -1.0, "floor(-0.5) == -1.0"); + check(floor(half) == 0.0 && !signbit(floor(half)), "floor(+0.5) == +0.0"); + check(floorl(neg_half) == -1.0L, "floorl(-0.5L) == -1.0L"); + check(floorf(-0.5f) == -1.0f, "floorf(-0.5f) == -1.0f"); + check(floor(two) == 2.0, "floor(2.0) == 2.0"); + check(floorf(twof) == 2.0f, "floorf(2.0f) == 2.0f"); + check(floor(-two) == -2.0, "floor(-2.0) == -2.0"); + check(floor(negzero) == 0.0 && signbit(floor(negzero)), "floor(-0.0) == -0.0"); + check(floor(sub) == 0.0 && !signbit(floor(sub)), "floor(5e-324) == +0.0"); + check(floor(-sub) == -1.0, "floor(-5e-324) == -1.0"); + check(floor(big_half) == big, "floor(2^51 + 0.5) == 2^51"); + check(floor(-big_half) == -big - 1.0, "floor(-(2^51 + 0.5)) == -2^51 - 1"); + check(floorl((long double)big + 0.5L) == (long double)big, "floorl(2^62 + 0.5L) == 2^62L"); + check(floorl(-((long double)big + 0.5L)) == -((long double)big) - 1.0L, + "floorl(-(2^62 + 0.5L)) == -2^62L - 1"); +} + +/* 4. ceil: toward +Inf; ceil(-0.5) must be a real negative zero. */ +static void +ceil_scenario(void) +{ + volatile double negzero = -0.0; + volatile double neg_half = -0.5; + volatile float pos_half = 0.5f; + volatile long double pos_half_l = 0.5L; + volatile double neg_third = -0.3; + volatile double big = 2251799813685248.0; + volatile double big_half = 2251799813685248.5; + + check(ceil(neg_half) == 0.0 && signbit(ceil(neg_half)), "ceil(-0.5) == -0.0 (signbit)"); + check(ceil(-1.5) == -1.0, "ceil(-1.5) == -1.0"); + check(ceilf(pos_half) == 1.0f, "ceilf(+0.5f) == 1.0f"); + check(ceill(pos_half_l) == 1.0L, "ceill(+0.5L) == 1.0L"); + check(ceil(neg_third) == 0.0 && signbit(ceil(neg_third)), "ceil(-0.3) == -0.0 (signbit)"); + check(ceil(negzero) == 0.0 && signbit(ceil(negzero)), "ceil(-0.0) == -0.0"); + check(ceilf(-0.5f) == 0.0f && signbit(ceilf(-0.5f)), "ceilf(-0.5f) == -0.0f (signbit)"); + check(ceill(-0.5L) == 0.0L && signbit(ceill(-0.5L)), "ceill(-0.5L) == -0.0L (signbit)"); + check(ceil(big_half) == big + 1.0, "ceil(2^51 + 0.5) == 2^51 + 1"); + check(ceil(-big_half) == -big, "ceil(-(2^51 + 0.5)) == -2^51"); +} + +/* 5. trunc: toward zero, with signed-zero results below 1. */ +static void +trunc_scenario(void) +{ + volatile double neg = -1.7; + volatile double pos = 1.7; + volatile double neg_tiny = -0.7; + volatile float posf = 1.7f; + volatile long double negl = -1.7L; + volatile double negzero = -0.0; + volatile double big = 2251799813685248.0; + volatile double big_half = 2251799813685248.5; + + check(trunc(neg) == -1.0, "trunc(-1.7) == -1.0"); + check(trunc(pos) == 1.0, "trunc(+1.7) == 1.0"); + check(truncf(posf) == 1.0f, "truncf(+1.7f) == 1.0f"); + check(truncl(negl) == -1.0L, "truncl(-1.7L) == -1.0L"); + check(trunc(neg_tiny) == 0.0 && signbit(trunc(neg_tiny)), "trunc(-0.7) == -0.0 (signbit)"); + check(trunc(0.7) == 0.0 && !signbit(trunc(0.7)), "trunc(+0.7) == +0.0"); + check(trunc(negzero) == 0.0 && signbit(trunc(negzero)), "trunc(-0.0) == -0.0"); + check(trunc(big_half) == big, "trunc(2^51 + 0.5) == 2^51"); + check(trunc(-big_half) == -big, "trunc(-(2^51 + 0.5)) == -2^51"); + check(truncl((long double)big + 0.5L) == (long double)big, "truncl(2^62 + 0.5L) == 2^62L"); +} + +/* 6. NaN and ±Inf passthrough for the three rounding families, plus the + * boundary value 4194304.5f (2^22 + 0.5, where float clearing drops a + * single bit). */ +static void +special_scenario(void) +{ + volatile double nan = NAN; + volatile double inf = HUGE_VAL; + volatile float nanf = NAN; + volatile long double inf_l = HUGE_VALL; + volatile float float_bound = 4194304.5f; + volatile double dbl_bound = 2251799813685248.5; + + check(isnan(floor(nan)), "floor(NaN) is NaN"); + check(isnan(ceil(nan)), "ceil(NaN) is NaN"); + check(isnan(trunc(nan)), "trunc(NaN) is NaN"); + check(isnan(floorf(nanf)), "floorf(NaN) is NaN"); + check(floor(inf) == inf, "floor(+Inf) == +Inf"); + check(floor(-inf) == -inf, "floor(-Inf) == -Inf"); + check(ceil(inf) == inf, "ceil(+Inf) == +Inf"); + check(ceil(-inf) == -inf, "ceil(-Inf) == -Inf"); + check(trunc(inf) == inf, "trunc(+Inf) == +Inf"); + check(trunc(-inf) == -inf, "trunc(-Inf) == -Inf"); + check(ceill(inf_l) == inf_l, "ceill(+Inf) == +Inf"); + check(truncl(-inf_l) == -inf_l, "truncl(-Inf) == -Inf"); + check(floorf(float_bound) == 4194304.0f, "floorf(2^22 + 0.5f) == 2^22f"); + check(trunc(float_bound) == 4194304.0, "trunc(2^22 + 0.5) == 2^22"); + check(ceilf(float_bound) == 4194305.0f, "ceilf(2^22 + 0.5f) == 2^22 + 1f"); + check(trunc(dbl_bound) == 2251799813685248.0, "trunc(2^51 + 0.5) == 2^51"); +} + +/* 7. Classification macros over the FP_* constants. */ +static void +classify_scenario(void) +{ + volatile double nan = NAN; + volatile double inf = HUGE_VAL; + volatile double zero = 0.0; + volatile double sub = 5e-324; + volatile double one = 1.0; + + check(fpclassify(nan) == FP_NAN, "fpclassify(NaN) == FP_NAN"); + check(fpclassify(inf) == FP_INFINITE, "fpclassify(+Inf) == FP_INFINITE"); + check(fpclassify(zero) == FP_ZERO, "fpclassify(0.0) == FP_ZERO"); + check(fpclassify(sub) == FP_SUBNORMAL, "fpclassify(5e-324) == FP_SUBNORMAL"); + check(fpclassify(one) == FP_NORMAL, "fpclassify(1.0) == FP_NORMAL"); + check(isnan(nan) && !isnan(one), "isnan distinguishes NaN"); + check(isinf(inf) && !isinf(one), "isinf distinguishes Inf"); + check(isfinite(one) && !isfinite(inf) && !isfinite(nan), "isfinite is false for Inf/NaN"); + check(isnormal(one) && !isnormal(zero) && !isnormal(sub) && !isnormal(nan), + "isnormal is false for zero/subnormal/NaN"); + check(signbit(-0.0) && !signbit(zero), "signbit distinguishes -0.0"); + check(fpclassify(1.0f) == FP_NORMAL, "fpclassify works on float"); + check(fpclassify(1.0L) == FP_NORMAL, "fpclassify works on long double"); +} + +int +main(void) +{ + fabs_scenario(); + copysign_scenario(); + floor_scenario(); + ceil_scenario(); + trunc_scenario(); + special_scenario(); + classify_scenario(); + + if (failures == 0) + { + say(1, "all math slice S1 tests passed\n"); + } + else + { + say(1, "FAILURES\n"); + } + return failures == 0 ? 0 : 1; +} diff --git a/tests/test_math_s2.c b/tests/test_math_s2.c new file mode 100644 index 0000000..e4eb8e2 --- /dev/null +++ b/tests/test_math_s2.c @@ -0,0 +1,357 @@ +/* + * vlibc — math.h slice S2 test: the rounding family (todo 39). + * + * Exercises every implemented family in all three precisions: + * + * 1. round/roundf/roundl — half away from zero: round(-0.5) == -1.0, + * round(2.5) == 3.0, signed zero preserved, subnormals collapse to + * ±0; + * 2. rint — round to nearest, ties to even, in the default (only + * reachable) rounding mode: rint(0.5) == 0.0, rint(-0.5) == -0.0 + * (a real negative zero), rint(2.5) == 2.0; + * 3. nearbyint — the same values as rint but through the guaranteed + * exception-free entry points; + * 4. lrint/lrintf/lrintl and llrint/llrintf/llrintl — the same nearest + * rounding returned as long / long long, over 0.0, integers, ties + * and the 2^40 scale (exact in both types); + * 5. ilogb/ilogbf/ilogbl — the unbiased exponent with the FP_ILOGB0 + * (zero) and FP_ILOGBNAN (Inf/NaN) sentinels, including the + * subnormal exponents -1074 (double), -149 (float) and -16445 + * (x87 long double); + * 6. logb/logbf/logbl — the same exponent returned as a floating + * value, with ±0 -> -Inf, ±Inf -> +Inf and NaN passing through. + * + * Expected values follow IEEE 754-2008 semantics by hand; they are NOT + * print-compared against the host libm at runtime (the host glibc on + * this machine returns FP_ILOGB0 for ilogb(NaN), deviating from C23, + * and the host TCB makes errno unreadable in this host-linked binary). + * Values that must survive exactly are fed through volatile objects so + * the compiler cannot constant-fold the call away; diagnostics go + * through raw SYS_write and no host headers are included. + */ + +#include + +#include "../src/internal/syscall.h" + +/* Write a NUL-terminated string to fd via the raw syscall layer. The + * optimize attribute keeps GCC from lowering the length loop into a + * strlen call, which would leave a vlibc-owned symbol undefined in this + * host-linked standalone binary (house idiom, see src/string). */ +static __attribute__((optimize("no-tree-loop-distribute-patterns"))) void +say(int fd, const char *s) +{ + long n = 0; + + while (s[n] != '\0') + { + n++; + } + __syscall3(SYS_write, fd, (long)s, n); +} + +static int failures; + +static void +check(int ok, const char *msg) +{ + if (ok) + { + say(1, "ok "); + } + else + { + say(1, "FAIL "); + failures++; + } + say(1, msg); + say(1, "\n"); +} + +/* 1. round: halfway cases away from zero, with signed zero and the + * subnormal collapse. */ +static void +round_scenario(void) +{ + volatile double negzero = -0.0; + volatile double sub = 0x1p-1074; /* smallest subnormal double */ + volatile double big = 1e300; + volatile long double ldsub = 0x1p-16445L; /* smallest subnormal ld */ + volatile float fsub = 0x1p-149f; /* smallest subnormal float */ + + check(round(0.0) == 0.0 && !signbit(round(0.0)), "round(+0.0) == +0.0"); + check(round(negzero) == 0.0 && signbit(round(negzero)), "round(-0.0) == -0.0"); + check(round(0.5) == 1.0, "round(0.5) == 1.0"); + check(round(-0.5) == -1.0, "round(-0.5) == -1.0"); + check(round(1.5) == 2.0, "round(1.5) == 2.0"); + check(round(-1.5) == -2.0, "round(-1.5) == -2.0"); + check(round(2.5) == 3.0, "round(2.5) == 3.0"); + check(round(-2.5) == -3.0, "round(-2.5) == -3.0"); + check(round(2.4) == 2.0, "round(2.4) == 2.0"); + check(round(2.6) == 3.0, "round(2.6) == 3.0"); + check(round(-2.4) == -2.0, "round(-2.4) == -2.0"); + check(round(-2.6) == -3.0, "round(-2.6) == -3.0"); + check(round(sub) == 0.0 && !signbit(round(sub)), "round(5e-324) == +0.0"); + check(round(-sub) == 0.0 && signbit(round(-sub)), "round(-5e-324) == -0.0"); + check(round(big) == big, "round(1e300) == 1e300"); + check(roundf(0.5f) == 1.0f, "roundf(0.5f) == 1.0f"); + check(roundf(-0.5f) == -1.0f, "roundf(-0.5f) == -1.0f"); + check(roundf(1.5f) == 2.0f, "roundf(1.5f) == 2.0f"); + check(roundf(-2.5f) == -3.0f, "roundf(-2.5f) == -3.0f"); + check(roundf(fsub) == 0.0f && !signbit(roundf(fsub)), "roundf(1e-45f) == +0.0f"); + check(roundf(-fsub) == 0.0f && signbit(roundf(-fsub)), "roundf(-1e-45f) == -0.0f"); + check(roundf(4194304.5f) == 4194305.0f, "roundf(2^22 + 0.5f) == 2^22 + 1f"); + check(roundl(0.5L) == 1.0L, "roundl(0.5L) == 1.0L"); + check(roundl(-0.5L) == -1.0L, "roundl(-0.5L) == -1.0L"); + check(roundl(2.5L) == 3.0L, "roundl(2.5L) == 3.0L"); + check(roundl(-2.5L) == -3.0L, "roundl(-2.5L) == -3.0L"); + check(roundl(ldsub) == 0.0L && !signbit(roundl(ldsub)), "roundl(min subnormal) == +0.0L"); + check(roundl(-ldsub) == 0.0L && signbit(roundl(-ldsub)), "roundl(-min subnormal) == -0.0L"); +} + +/* 2. rint: round to nearest, ties to even. */ +static void +rint_scenario(void) +{ + volatile double negzero = -0.0; + volatile double sub = 0x1p-1074; + volatile double big = 1e300; + volatile double half_int = 0x1p52 + 0.5; /* 2^52 + 0.5, tie to even */ + volatile long double ldsub = 0x1p-16445L; + volatile float fsub = 0x1p-149f; + + check(rint(0.0) == 0.0 && !signbit(rint(0.0)), "rint(+0.0) == +0.0"); + check(rint(negzero) == 0.0 && signbit(rint(negzero)), "rint(-0.0) == -0.0"); + check(rint(0.5) == 0.0 && !signbit(rint(0.5)), "rint(0.5) == +0.0 (tie to even)"); + check(rint(-0.5) == 0.0 && signbit(rint(-0.5)), "rint(-0.5) == -0.0 (tie to even)"); + check(rint(1.5) == 2.0, "rint(1.5) == 2.0"); + check(rint(-1.5) == -2.0, "rint(-1.5) == -2.0"); + check(rint(2.5) == 2.0, "rint(2.5) == 2.0 (tie to even)"); + check(rint(-2.5) == -2.0, "rint(-2.5) == -2.0 (tie to even)"); + check(rint(2.4) == 2.0, "rint(2.4) == 2.0"); + check(rint(2.6) == 3.0, "rint(2.6) == 3.0"); + check(rint(-2.6) == -3.0, "rint(-2.6) == -3.0"); + check(rint(sub) == 0.0 && !signbit(rint(sub)), "rint(5e-324) == +0.0"); + check(rint(-sub) == 0.0 && signbit(rint(-sub)), "rint(-5e-324) == -0.0"); + check(rint(big) == big, "rint(1e300) == 1e300"); + check(rint(half_int) == 0x1p52, "rint(2^52 + 0.5) == 2^52 (tie to even)"); + check(rint(-half_int) == -(0x1p52), "rint(-(2^52 + 0.5)) == -2^52"); + check(rintf(0.5f) == 0.0f && !signbit(rintf(0.5f)), "rintf(0.5f) == +0.0f"); + check(rintf(-0.5f) == 0.0f && signbit(rintf(-0.5f)), "rintf(-0.5f) == -0.0f"); + check(rintf(1.5f) == 2.0f, "rintf(1.5f) == 2.0f"); + check(rintf(2.5f) == 2.0f, "rintf(2.5f) == 2.0f (tie to even)"); + check(rintf(-2.5f) == -2.0f, "rintf(-2.5f) == -2.0f"); + check(rintf(4194304.5f) == 4194304.0f, "rintf(2^22 + 0.5f) == 2^22f (tie to even)"); + check(rintf(fsub) == 0.0f && !signbit(rintf(fsub)), "rintf(1e-45f) == +0.0f"); + check(rintl(0.5L) == 0.0L && !signbit(rintl(0.5L)), "rintl(0.5L) == +0.0L"); + check(rintl(-0.5L) == 0.0L && signbit(rintl(-0.5L)), "rintl(-0.5L) == -0.0L"); + check(rintl(1.5L) == 2.0L, "rintl(1.5L) == 2.0L"); + check(rintl(2.5L) == 2.0L, "rintl(2.5L) == 2.0L (tie to even)"); + check(rintl(-2.5L) == -2.0L, "rintl(-2.5L) == -2.0L"); + check(rintl(ldsub) == 0.0L && !signbit(rintl(ldsub)), "rintl(min subnormal) == +0.0L"); +} + +/* 3. nearbyint: same nearest-even rounding, exception-free entry. */ +static void +nearbyint_scenario(void) +{ + volatile double negzero = -0.0; + volatile double sub = 0x1p-1074; + volatile long double ldsub = 0x1p-16445L; + volatile float fsub = 0x1p-149f; + + check(nearbyint(0.0) == 0.0 && !signbit(nearbyint(0.0)), "nearbyint(+0.0) == +0.0"); + check(nearbyint(negzero) == 0.0 && signbit(nearbyint(negzero)), "nearbyint(-0.0) == -0.0"); + check(nearbyint(0.5) == 0.0 && !signbit(nearbyint(0.5)), "nearbyint(0.5) == +0.0"); + check(nearbyint(-0.5) == 0.0 && signbit(nearbyint(-0.5)), "nearbyint(-0.5) == -0.0"); + check(nearbyint(1.5) == 2.0, "nearbyint(1.5) == 2.0"); + check(nearbyint(2.5) == 2.0, "nearbyint(2.5) == 2.0 (tie to even)"); + check(nearbyint(-2.5) == -2.0, "nearbyint(-2.5) == -2.0"); + check(nearbyint(2.4) == 2.0, "nearbyint(2.4) == 2.0"); + check(nearbyint(2.6) == 3.0, "nearbyint(2.6) == 3.0"); + check(nearbyint(-2.6) == -3.0, "nearbyint(-2.6) == -3.0"); + check(nearbyint(sub) == 0.0 && !signbit(nearbyint(sub)), "nearbyint(5e-324) == +0.0"); + check(nearbyint(-sub) == 0.0 && signbit(nearbyint(-sub)), "nearbyint(-5e-324) == -0.0"); + check(nearbyintf(0.5f) == 0.0f, "nearbyintf(0.5f) == +0.0f"); + check(nearbyintf(2.5f) == 2.0f, "nearbyintf(2.5f) == 2.0f"); + check(nearbyintf(-0.5f) == 0.0f && signbit(nearbyintf(-0.5f)), "nearbyintf(-0.5f) == -0.0f"); + check(nearbyintf(fsub) == 0.0f, "nearbyintf(1e-45f) == +0.0f"); + check(nearbyintl(0.5L) == 0.0L, "nearbyintl(0.5L) == +0.0L"); + check(nearbyintl(-0.5L) == 0.0L && signbit(nearbyintl(-0.5L)), "nearbyintl(-0.5L) == -0.0L"); + check(nearbyintl(2.5L) == 2.0L, "nearbyintl(2.5L) == 2.0L"); + check(nearbyintl(-2.5L) == -2.0L, "nearbyintl(-2.5L) == -2.0L"); + check(nearbyintl(ldsub) == 0.0L, "nearbyintl(min subnormal) == +0.0L"); +} + +/* 4. lrint/llrint: nearest rounding narrowed to long / long long. */ +static void +lrint_scenario(void) +{ + volatile double two40 = 0x1p40; /* 1099511627776.0, exact in long */ + + check(lrint(0.0) == 0L, "lrint(0.0) == 0"); + check(lrint(1.0) == 1L, "lrint(1.0) == 1"); + check(lrint(1.5) == 2L, "lrint(1.5) == 2"); + check(lrint(2.5) == 2L, "lrint(2.5) == 2 (tie to even)"); + check(lrint(-1.5) == -2L, "lrint(-1.5) == -2"); + check(lrint(-2.5) == -2L, "lrint(-2.5) == -2"); + check(lrint(0.5) == 0L, "lrint(0.5) == 0"); + check(lrint(-0.5) == 0L, "lrint(-0.5) == 0"); + check(lrint(two40) == 1099511627776L, "lrint(2^40) == 2^40"); + check(lrintf(1.5f) == 2L, "lrintf(1.5f) == 2"); + check(lrintf(2.5f) == 2L, "lrintf(2.5f) == 2"); + check(lrintf(-1.5f) == -2L, "lrintf(-1.5f) == -2"); + check(lrintf(0.0f) == 0L, "lrintf(0.0f) == 0"); + check(lrintl(1.5L) == 2L, "lrintl(1.5L) == 2"); + check(lrintl(2.5L) == 2L, "lrintl(2.5L) == 2"); + check(lrintl(-2.5L) == -2L, "lrintl(-2.5L) == -2"); + check(lrintl(0.0L) == 0L, "lrintl(0.0L) == 0"); + check(lrintl(-1.5L) == -2L, "lrintl(-1.5L) == -2"); +} + +static void +llrint_scenario(void) +{ + volatile double two40 = 0x1p40; + + check(llrint(0.0) == 0LL, "llrint(0.0) == 0"); + check(llrint(1.0) == 1LL, "llrint(1.0) == 1"); + check(llrint(1.5) == 2LL, "llrint(1.5) == 2"); + check(llrint(2.5) == 2LL, "llrint(2.5) == 2 (tie to even)"); + check(llrint(-1.5) == -2LL, "llrint(-1.5) == -2"); + check(llrint(-2.5) == -2LL, "llrint(-2.5) == -2"); + check(llrint(two40) == 1099511627776LL, "llrint(2^40) == 2^40"); + check(llrintf(1.5f) == 2LL, "llrintf(1.5f) == 2"); + check(llrintf(-2.5f) == -2LL, "llrintf(-2.5f) == -2"); + check(llrintf(0.0f) == 0LL, "llrintf(0.0f) == 0"); + check(llrintl(2.5L) == 2LL, "llrintl(2.5L) == 2"); + check(llrintl(-1.5L) == -2LL, "llrintl(-1.5L) == -2"); + check(llrintl(1.5L) == 2LL, "llrintl(1.5L) == 2"); +} + +/* 5. ilogb: floor(log2 |x|) with the zero and Inf/NaN sentinels. */ +static void +ilogb_scenario(void) +{ + volatile double nan = NAN; + volatile double inf = HUGE_VAL; + volatile double negzero = -0.0; + volatile float nanf = NAN; + volatile long double nana_l = NAN; + volatile long double ldsub = 0x1p-16445L; + volatile float fsub = 0x1p-149f; + + check(ilogb(0.0) == FP_ILOGB0, "ilogb(+0.0) == FP_ILOGB0"); + check(ilogb(negzero) == FP_ILOGB0, "ilogb(-0.0) == FP_ILOGB0"); + check(ilogb(nan) == FP_ILOGBNAN, "ilogb(NaN) == FP_ILOGBNAN"); + check(ilogb(inf) == FP_ILOGBNAN, "ilogb(+Inf) == FP_ILOGBNAN"); + check(ilogb(-inf) == FP_ILOGBNAN, "ilogb(-Inf) == FP_ILOGBNAN"); + check(ilogb(1.0) == 0, "ilogb(1.0) == 0"); + check(ilogb(2.0) == 1, "ilogb(2.0) == 1"); + check(ilogb(8.0) == 3, "ilogb(8.0) == 3"); + check(ilogb(-8.0) == 3, "ilogb(-8.0) == 3"); + check(ilogb(0.5) == -1, "ilogb(0.5) == -1"); + check(ilogb(0x1p-1074) == -1074, "ilogb(5e-324) == -1074"); + check(ilogb(0x1p52) == 52, "ilogb(2^52) == 52"); + check(ilogb(0x1p53) == 53, "ilogb(2^53) == 53"); + check(ilogbf(1.0f) == 0, "ilogbf(1.0f) == 0"); + check(ilogbf(8.0f) == 3, "ilogbf(8.0f) == 3"); + check(ilogbf(0.5f) == -1, "ilogbf(0.5f) == -1"); + check(ilogbf(0.0f) == FP_ILOGB0, "ilogbf(0.0f) == FP_ILOGB0"); + check(ilogbf(nanf) == FP_ILOGBNAN, "ilogbf(NaN) == FP_ILOGBNAN"); + check(ilogbf(fsub) == -149, "ilogbf(min subnormal) == -149"); + check(ilogbl(1.0L) == 0, "ilogbl(1.0L) == 0"); + check(ilogbl(8.0L) == 3, "ilogbl(8.0L) == 3"); + check(ilogbl(0.5L) == -1, "ilogbl(0.5L) == -1"); + check(ilogbl(0.0L) == FP_ILOGB0, "ilogbl(0.0L) == FP_ILOGB0"); + check(ilogbl(nana_l) == FP_ILOGBNAN, "ilogbl(NaN) == FP_ILOGBNAN"); + check(ilogbl(ldsub) == -16445, "ilogbl(min subnormal) == -16445"); +} + +/* 6. logb: the same exponent as a floating value. */ +static void +logb_scenario(void) +{ + volatile double nan = NAN; + volatile double inf = HUGE_VAL; + volatile double negzero = -0.0; + volatile float nanf = NAN; + volatile long double nana_l = NAN; + volatile long double ldsub = 0x1p-16445L; + volatile float fsub = 0x1p-149f; + + check(logb(8.0) == 3.0, "logb(8.0) == 3.0"); + check(logb(-8.0) == 3.0, "logb(-8.0) == 3.0"); + check(logb(0.5) == -1.0, "logb(0.5) == -1.0"); + check(logb(0x1p-1074) == -1074.0, "logb(5e-324) == -1074.0"); + check(logb(inf) == inf, "logb(+Inf) == +Inf"); + check(logb(-inf) == inf, "logb(-Inf) == +Inf"); + check(isnan(logb(nan)), "logb(NaN) is NaN"); + check(logb(0.0) == -HUGE_VAL, "logb(+0.0) == -Inf"); + check(logb(negzero) == -HUGE_VAL, "logb(-0.0) == -Inf"); + check(logbf(8.0f) == 3.0f, "logbf(8.0f) == 3.0f"); + check(logbf(0.5f) == -1.0f, "logbf(0.5f) == -1.0f"); + check(logbf(fsub) == -149.0f, "logbf(min subnormal) == -149.0f"); + check(logbf(0.0f) == -HUGE_VALF, "logbf(0.0f) == -Inf"); + check(logbf(1.0f) == 0.0f, "logbf(1.0f) == 0.0f"); + check(isnan(logbf(nanf)), "logbf(NaN) is NaN"); + check(logbl(0.5L) == -1.0L, "logbl(0.5L) == -1.0L"); + check(logbl(ldsub) == -16445.0L, "logbl(min subnormal) == -16445.0L"); + check(logbl(0.0L) == -HUGE_VALL, "logbl(0.0L) == -Inf"); + check(logbl(1.0L) == 0.0L, "logbl(1.0L) == 0.0L"); + check(logbl(8.0L) == 3.0L, "logbl(8.0L) == 3.0L"); + check(isnan(logbl(nana_l)), "logbl(NaN) is NaN"); +} + +/* 7. NaN and ±Inf passthrough for the three rounding families. */ +static void +special_scenario(void) +{ + volatile double nan = NAN; + volatile double inf = HUGE_VAL; + volatile float nanf = NAN; + volatile long double nana_l = NAN; + volatile long double inf_l = HUGE_VALL; + + check(isnan(round(nan)), "round(NaN) is NaN"); + check(isnan(rint(nan)), "rint(NaN) is NaN"); + check(isnan(nearbyint(nan)), "nearbyint(NaN) is NaN"); + check(round(inf) == inf, "round(+Inf) == +Inf"); + check(round(-inf) == -inf, "round(-Inf) == -Inf"); + check(rint(inf) == inf, "rint(+Inf) == +Inf"); + check(rint(-inf) == -inf, "rint(-Inf) == -Inf"); + check(nearbyint(inf) == inf, "nearbyint(+Inf) == +Inf"); + check(nearbyint(-inf) == -inf, "nearbyint(-Inf) == -Inf"); + check(isnan(roundf(nanf)), "roundf(NaN) is NaN"); + check(isnan(rintf(nanf)), "rintf(NaN) is NaN"); + check(isnan(nearbyintf(nanf)), "nearbyintf(NaN) is NaN"); + check(isnan(roundl(nana_l)), "roundl(NaN) is NaN"); + check(isnan(rintl(nana_l)), "rintl(NaN) is NaN"); + check(isnan(nearbyintl(nana_l)), "nearbyintl(NaN) is NaN"); + check(roundl(inf_l) == inf_l, "roundl(+Inf) == +Inf"); + check(rintl(-inf_l) == -inf_l, "rintl(-Inf) == -Inf"); + check(nearbyintl(inf_l) == inf_l, "nearbyintl(+Inf) == +Inf"); +} + +int +main(void) +{ + round_scenario(); + rint_scenario(); + nearbyint_scenario(); + lrint_scenario(); + llrint_scenario(); + ilogb_scenario(); + logb_scenario(); + special_scenario(); + + if (failures == 0) + { + say(1, "all math slice S2 tests passed\n"); + } + else + { + say(1, "FAILURES\n"); + } + return failures == 0 ? 0 : 1; +}