math.hpp File Reference

Floating-point functions that give bit-identical results on every compiler, standard library and platform. More...

#include <array>
#include <cfloat>
#include <cmath>
#include <cstddef>
#include <limits>
#include <numbers>

Go to the source code of this file.

Functions
double	dasmig::core::math::fenced (double x) noexcept
	Return x unchanged, hiding its value from the optimizer.
double	dasmig::core::math::exp (double x) noexcept
	Natural exponential, bit-identical on every platform.
double	dasmig::core::math::log (double x) noexcept
	Natural logarithm, bit-identical on every platform.

Detailed Description

Floating-point functions that give bit-identical results on every compiler, standard library and platform.

Author

Diego Dasso Migotto (diegomigotto at hotmail dot com)

std::exp and std::log are not required to be correctly rounded, and glibc, the MSVC runtime, Apple's libm and Emscripten's musl differ in the last bit for some inputs. The functions here use only operations IEEE 754 defines exactly (+, -, *, /, sqrt, floor, frexp, ldexp), so the same input gives the same bits everywhere. They are accurate to a few ulp, which is what sampling needs; they are not replacements for <cmath> in general.

Exactness also needs every operation rounded on its own. Compilers may fuse a * b + c into one FMA instruction (GCC in its default GNU mode, clang with -ffp-contract=fast, ARM64 targets), which rounds once instead of twice. Pragmas do not prevent this reliably (clang ignores them under -ffp-contract=fast), so every product that feeds an addition in the core passes through fenced(), which the compiler cannot see through. Fast-math modes, which also reassociate, are refused at compile time.

Definition in file math.hpp.

Function Documentation

exp()

double dasmig::core::math::exp ( double  x )

inlinenoexcept

Natural exponential, bit-identical on every platform.

Reduces x = k ln(2) + r with |r| <= ln(2)/2, evaluates a degree-13 Taylor polynomial in r by Horner's rule and scales by 2^k.

Parameters

x

Exponent.

Returns

e^x; +inf above about 709.78, 0 below about -745.13, NaN for NaN.

Definition at line 128 of file math.hpp.

fenced()

double dasmig::core::math::fenced ( double  x )

inlinenoexcept

Return x unchanged, hiding its value from the optimizer.

Wrap a product in this before adding to it, so the compiler cannot fuse the two into a fused multiply-add: fenced(a * b) + c. Costs nothing on x86-64 and ARM64 (an empty asm statement); WebAssembly has no scalar FMA, and MSVC (VS2022, /fp:precise) does not contract, so neither needs it.

Parameters

x	A value, typically a product.

Returns

x, rounded to double.

Definition at line 60 of file math.hpp.

log()

double dasmig::core::math::log ( double  x )

inlinenoexcept

Natural logarithm, bit-identical on every platform.

Splits x = m 2^e with m in [sqrt(1/2), sqrt(2)), evaluates log(m) = 2 atanh((m-1)/(m+1)) as an odd series and adds e ln(2).

Parameters

x

Argument.

Returns

ln(x); -inf for 0, +inf for +inf, NaN for negative or NaN.

Definition at line 154 of file math.hpp.