vlibc/include/math.h

#ifndef VLIBC_MATH_H
#define VLIBC_MATH_H

/*
 * vlibc — <math.h>.
 *
 * 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 <complex.h> 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 <vlibc/features.h>

#include <limits.h>
#include <stddef.h>

#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
 * <limits.h> 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 */