diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md index 9ec86ef4..57fdb035 100644 --- a/.github/copilot-instructions.md +++ b/.github/copilot-instructions.md @@ -261,3 +261,5 @@ See existing operations like `Mul`, `Add`, `SaturateCast` in `include/fused_kern 1. **Windows Ninja + NVCC path**: After CMake configure on Windows with Ninja, `/CMakeFiles/rules.ninja` may contain an incorrect path to `nvcc.exe`. The CI workflow patches this with PowerShell `Set-Content`. If you hit this locally, check that `CUDACXX` env var is set before invoking CMake and verify the generated `rules.ninja`. +2. **CUDA Compilation Availability**: CUDA compilation is supported and should not be declared blocked based on an unsuccessful ad hoc nvcc invocation; use the repository's documented VS Developer Shell/CMake build configuration and correct architecture/toolchain settings before concluding CUDA is unavailable. + diff --git a/include/fused_kernel/core/constexpr_libs/constexpr_cmath.h b/include/fused_kernel/core/constexpr_libs/constexpr_cmath.h index 95e538e6..766fb984 100644 --- a/include/fused_kernel/core/constexpr_libs/constexpr_cmath.h +++ b/include/fused_kernel/core/constexpr_libs/constexpr_cmath.h @@ -27,6 +27,7 @@ #ifdef __CUDACC__ #include #include +#include // Conditionally include the algorithm header if compiling on CUDA 13.3+ #if __has_include() @@ -43,6 +44,7 @@ using cuda::std::cmp_less; using cuda::std::cmp_less_equal; using cuda::std::cmp_not_equal; using cuda::std::is_constant_evaluated; +using cuda::std::numeric_limits; // If the header exists, alias the cuda::std versions #if __has_include() @@ -78,23 +80,24 @@ using std::cmp_not_equal; using std::max; using std::min; using std::is_constant_evaluated; +using std::numeric_limits; } // namespace base } // namespace cxp #endif namespace cxp { template - constexpr T minValue = std::numeric_limits::lowest(); + constexpr T minValue = base::numeric_limits::lowest(); template - constexpr T maxValue = std::numeric_limits::max(); + constexpr T maxValue = base::numeric_limits::max(); template - constexpr T smallestPositiveValue = std::is_floating_point_v ? std::numeric_limits::min() : static_cast(1); + constexpr T smallestPositiveValue = std::is_floating_point_v ? base::numeric_limits::min() : static_cast(1); #define CXP_F_FUNC \ template \ - FK_HOST_DEVICE_FUSE auto f(const Types... vals) { \ + FK_HOST_DEVICE_FUSE decltype(auto) f(const Types... vals) { \ return Exec::exec(vals...); \ } @@ -666,7 +669,7 @@ namespace cxp { { if (base::is_constant_evaluated()) { // Replace union with standard C++20 bit_cast - uint ui = cxp::bit_cast(x); + uint ui = bit_cast(x); // Extract the 8-bit exponent int32_t e = (ui >> 23) & 0xFF; @@ -678,7 +681,7 @@ namespace cxp { return x; // Subnormal normalization: multiply by 2^24 to push into normal range - ui = cxp::bit_cast(x * 16777216.0f); // Re-cast after math + ui = bit_cast(x * 16777216.0f); // Re-cast after math e = ((ui >> 23) & 0xFF) - 24; } else if (e == 0xFF) { // Infinity or NaN @@ -692,7 +695,7 @@ namespace cxp { if (e > 254) { // Force the exponent to the infinity marker while keeping the sign bit ui = (ui & 0x80000000) | 0x7F800000; - return cxp::bit_cast(ui); + return bit_cast(ui); } // 4. Check for Underflow @@ -700,12 +703,12 @@ namespace cxp { // Total underflow to zero if (e <= -24) { ui &= 0x80000000; // Preserve sign bit, zero everything else - return cxp::bit_cast(ui); + return bit_cast(ui); } // Partial underflow to subnormal ui = (ui & 0x807FFFFF) | ((e + 24) << 23); - return cxp::bit_cast(ui) * 5.960464477539063e-8f; + return bit_cast(ui) * 5.960464477539063e-8f; } // 5. Standard Reconstruction @@ -721,76 +724,101 @@ namespace cxp { CXP_F_FUNC }; - struct expf { + struct exp { struct BaseFunc { using InstanceType = fk::UnaryType; - FK_HOST_DEVICE_FUSE auto exec(const float x) { + template + FK_HOST_DEVICE_FUSE ST exec(const ST x) { if (base::is_constant_evaluated()) { - // 1. Handle edge cases FIRST to protect constexpr evaluation. - // The finite cutoffs sit safely outside the domain where the result is still - // representable, so they only fire where the answer is exactly inf or zero. - // Narrower bounds would misreport values near the overflow boundary, whose - // true result is still finite. - if (isnan::BaseFunc::exec(x)) + return c_exp(x); + } else { + return std::exp(x); + } + } + + private: + template + FK_HOST_DEVICE_FUSE ST c_exp(const ST x) { + if constexpr (std::is_same_v) { + if (x != x) return x; - if (x > 89.0f) - return cxp::bit_cast(0x7F800000); // 0x7F800000 is +INFINITY if (x < -104.0f) - return 0.0f; // Hard underflow, below half the smallest subnormal - - // The whole reduction is carried out in double and rounded to float exactly once - // at the end. float has 24 mantissa bits against double's 53, so the ~29 bits of - // headroom make the single final rounding correct, which is what makes this path - // agree with std::exp on a float input. - const double xd = static_cast(x); - - // 2. Range Reduction: x = k * ln(2) + r, with |r| <= ln(2)/2 - const double INV_LN2 = 1.4426950408889634; - const double kd = xd * INV_LN2; - // Round half away from zero without leaving constant evaluation - const int32_t k = static_cast(kd >= 0.0 ? kd + 0.5 : kd - 0.5); - const double k_d = static_cast(k); - - // ln(2) split so that k * LN2_HI is exact for every k we can reach here - const double LN2_HI = 6.93147180369123816490e-01; - const double LN2_LO = 1.90821492927058770002e-10; - const double r = (xd - k_d * LN2_HI) - k_d * LN2_LO; - - // 3. Degree-13 Taylor series in double. The truncation error over |r| <= 0.3466 - // is below 0.02 double ULP, far under half a float ULP. - const double poly = - 1.0 + - r * (1.0 + - r * (1.0 / 2.0 + - r * (1.0 / 6.0 + - r * (1.0 / 24.0 + - r * (1.0 / 120.0 + - r * (1.0 / 720.0 + - r * (1.0 / 5040.0 + - r * (1.0 / 40320.0 + - r * (1.0 / 362880.0 + - r * (1.0 / 3628800.0 + - r * (1.0 / 39916800.0 + - r * (1.0 / 479001600.0 + - r * (1.0 / - 6227020800.0))))))))))))); - - // 4. Reconstruction: e^x = e^r * 2^k. k stays well inside the double exponent - // range, so this scaling is exact and cannot underflow before the final rounding. - const double twoK = cxp::bit_cast(static_cast(k + 1023) << 52); - const double result = poly * twoK; - - // Anything at or above the float rounding midpoint towards 2^128 becomes infinity. - // Guarding here keeps the narrowing conversion below in range. - if (result >= 3.4028235677973366e38) { - return cxp::bit_cast(0x7F800000); + return 0.0f; + // 88.7228394f (0x42B17218) is the first input whose exponential overflows. + // The last finite input is its predecessor, 88.7228088f (0x42B17214). + if (x >= bit_cast(0x42B17218u)) + return base::numeric_limits::infinity(); + + float y = x * 1.44269504089f; + int k = static_cast((y >= 0.0f) ? (y + 0.5f) : (y - 0.5f)); + + // Split LN2 for float + float r = (x - static_cast(k) * 0.693145751953125f) - + static_cast(k) * 1.428606765330187e-06f; + + float r2 = r * r; + float poly = + 1.0f + r + + r2 * (0.5f + r * (0.16666667163f + + r * (0.04166666790f + + r * (0.00833333333f + r * (0.00138888889f + r * 0.00019841269f))))); + + if (k <= -127) { + // 2^k is not representable as a normal float: delegate to ldexpf, + // which performs correctly rounded gradual underflow to subnormals. + return ldexpf::BaseFunc::exec(poly, k); } - // Single, correctly rounded narrowing. Underflow to (sub)normal or zero is - // well defined and also correctly rounded. - return static_cast(result); + if (k >= 128) { + // 2^k is not representable: scale by 2^127 and then by 2. + // Both multiplications are exact power-of-two scalings, so no + // extra rounding occurs and overflow saturates to +inf naturally. + return poly * bit_cast(static_cast(254) << 23) * 2.0f; + } + + uint pow2_bits = static_cast(k + 127) << 23; + return poly * bit_cast(pow2_bits); + } else { + if (x != x) + return x; + if (x < -709.0) + return 0.0; + if (x > 709.0) + return base::numeric_limits::infinity(); + + double y = x * 1.44269504088896340736; + long long k = static_cast((y >= 0.0) ? (y + 0.5) : (y - 0.5)); + + double r = (x - static_cast(k) * 0.69314718036912381649) - + static_cast(k) * 1.90821492927058770002e-10; + + double r2 = r * r; + double poly = 1.0 + r + r2 * (1.0 / 2.0 + + r * (1.0 / 6.0 + + r * (1.0 / 24.0 + + r * (1.0 / 120.0 + + r * (1.0 / 720.0 + + r * (1.0 / 5040.0 + + r * (1.0 / 40320.0 + + r * (1.0 / 362880.0 + + r * (1.0 / 3628800.0 + + r * (1.0 / 39916800.0 + + r * 1.0 / 479001600.0)))))))))); + + ulonglong pow2_bits = static_cast(k + 1023) << 52; + return poly * bit_cast(pow2_bits); + } + } + }; + CXP_F_FUNC + }; + + struct expf { + struct BaseFunc { + using InstanceType = fk::UnaryType; + FK_HOST_DEVICE_FUSE float exec(const float x) { + if (base::is_constant_evaluated()) { + return exp::BaseFunc::exec(x); } else { - // The constexpr body was verified to be bit - // identical to std::exp for every finite float in [-104, 89]. return std::exp(x); } } @@ -812,6 +840,216 @@ namespace cxp { } }; + // --------------------------------------------------------- + // Accurate constexpr natural logarithm (c_logf) + // --------------------------------------------------------- + struct log { + struct BaseFunc { + using InstanceType = fk::UnaryType; + template + FK_HOST_DEVICE_FUSE ST exec(const ST x) { + if (base::is_constant_evaluated()) { + return c_log(x); + } else { + return std::log(x); + } + } + + private: + template + FK_HOST_DEVICE_FUSE ST c_log(const ST x) { + if constexpr (std::is_same_v) { + if (x != x) + return x; + if (x < 0.0f) + return base::numeric_limits::quiet_NaN(); + if (x == 0.0f) + return -base::numeric_limits::infinity(); + if (x == base::numeric_limits::infinity()) + return x; + + uint bits = bit_cast(x); + int e = static_cast((bits >> 23) & 0xFF) - 127; + + float m = bit_cast((bits & 0x007FFFFF) | 0x3F800000); + if (m > 1.41421356237f) { + m /= 2.0f; + e += 1; + } + + float z = (m - 1.0f) / (m + 1.0f); + float z2 = z * z; + + // 32-bit Remez Minimax Coefficients + float p = + z * (2.0f + z2 * (0.66666662693f + + z2 * (0.39999997616f + + z2 * (0.28571429849f + + z2 * (0.22222198546f + + z2 * (0.18183572590f + + z2 * (0.15313838422f + + z2 * 0.14798198640f))))))); + + return p + static_cast(e) * 0.6931471805f; + } else { + if (x != x) + return x; + if (x < 0.0) + return base::numeric_limits::quiet_NaN(); + if (x == 0.0) + return -base::numeric_limits::infinity(); + if (x == base::numeric_limits::infinity()) + return x; + + ulonglong bits = bit_cast(x); + int e = static_cast(bits >> 52) - 1023; + double m = bit_cast((bits & 0x000FFFFFFFFFFFFFull) | 0x3FF0000000000000ull); + + if (m > 1.4142135623730950488) { + m /= 2.0; + e += 1; + } + + double z = (m - 1.0) / (m + 1.0); + double z2 = z * z; + + // 64-bit Remez Minimax Coefficients + double p = z * (2.0 + z2 * (0.6666666666666735130 + + z2 * (0.3999999999940941908 + + z2 * (0.2857142874366239149 + + z2 * (0.2222219843214978396 + + z2 * (0.1818357216161805012 + + z2 * (0.1531383769920937332 + + z2 * 0.1479819860511658591))))))); + + return p + static_cast(e) * 0.693147180559945309417; + } + } + }; + CXP_F_FUNC + }; + + struct logf { + struct BaseFunc { + using InstanceType = fk::UnaryType; + FK_HOST_DEVICE_FUSE float exec(const float x) { + if (base::is_constant_evaluated()) { + return log::BaseFunc::exec(x); + } else { + return std::log(x); + } + } + }; + CXP_F_FUNC + }; + + struct pow { + struct BaseFunc { + using InstanceType = fk::BinaryType; + template + FK_HOST_DEVICE_FUSE ST exec(const ST base, const ST exponent) { + if (base::is_constant_evaluated()) { + return c_pow(base, exponent); + } else { + return std::pow(base, exponent); + } + } + + private: + FK_HOST_DEVICE_FUSE void split_float(float a, float &hi, float &lo) { + constexpr float split_factor = 4097.0f; // (1 << 12) + 1.0f + float c = split_factor * a; + float abig = c - a; + hi = c - abig; + lo = a - hi; + } + + FK_HOST_DEVICE_FUSE void split_double(double a, double &hi, double &lo) { + constexpr double split_factor = 134217729.0; // (1ull << 27) + 1.0 + double c = split_factor * a; + double abig = c - a; + hi = c - abig; + lo = a - hi; + } + + template + FK_HOST_DEVICE_FUSE ST c_pow(const ST base, const ST exponent) { + if constexpr (std::is_same_v) { + if (exponent == 0.0f || base == 1.0f) + return 1.0f; + if (base != base || exponent != exponent) + return base::numeric_limits::quiet_NaN(); + if (base == 0.0f) + return (exponent > 0.0f) ? 0.0f : base::numeric_limits::infinity(); + + if (base < 0.0f) { + int e_int = static_cast(exponent); + if (static_cast(e_int) == exponent) { + float res = exp::BaseFunc::exec(exponent * log::BaseFunc::exec(-base)); + return (e_int % 2 != 0) ? -res : res; + } + return base::numeric_limits::quiet_NaN(); + } + + float ln_x = log::BaseFunc::exec(base); + + // 24-bit mantissa double-split arithmetic (simulating 48-bit float precision) + float y_hi, y_lo, ln_x_hi, ln_x_lo; + split_float(exponent, y_hi, y_lo); + split_float(ln_x, ln_x_hi, ln_x_lo); + + float prod = exponent * ln_x; + float prod_err = ((y_hi * ln_x_hi - prod) + y_hi * ln_x_lo + y_lo * ln_x_hi) + y_lo * ln_x_lo; + + return exp::BaseFunc::exec(prod + prod_err); + } else { + if (exponent == 0.0 || base == 1.0) + return 1.0; + if (base != base || exponent != exponent) + return base::numeric_limits::quiet_NaN(); + if (base == 0.0) + return (exponent > 0.0) ? 0.0 : base::numeric_limits::infinity(); + + if (base < 0.0) { + long long e_int = static_cast(exponent); + if (static_cast(e_int) == exponent) { + double res = exp::BaseFunc::exec(exponent * log::BaseFunc::exec(-base)); + return (e_int % 2 != 0) ? -res : res; + } + return base::numeric_limits::quiet_NaN(); + } + + double ln_x = log::BaseFunc::exec(base); + + // 53-bit mantissa double-split arithmetic (simulating 106-bit double precision) + double y_hi, y_lo, ln_x_hi, ln_x_lo; + split_double(exponent, y_hi, y_lo); + split_double(ln_x, ln_x_hi, ln_x_lo); + + double prod = exponent * ln_x; + double prod_err = ((y_hi * ln_x_hi - prod) + y_hi * ln_x_lo + y_lo * ln_x_hi) + y_lo * ln_x_lo; + + return exp::BaseFunc::exec(prod + prod_err); + } + } + }; + CXP_F_FUNC + }; + + struct powf { + struct BaseFunc { + using InstanceType = fk::BinaryType; + FK_HOST_DEVICE_FUSE float exec(const float base, const float exponent) { + if (base::is_constant_evaluated()) { + return pow::BaseFunc::exec(base, exponent); + } else { + return std::pow(base, exponent); + } + } + }; + CXP_F_FUNC + }; + #undef CXP_F_FUNC } // namespace cxp diff --git a/utests/core/constexpr_libs/utest_constexpr_cmath.h b/utests/core/constexpr_libs/utest_constexpr_cmath.h index df4303fe..9c59a32b 100644 --- a/utests/core/constexpr_libs/utest_constexpr_cmath.h +++ b/utests/core/constexpr_libs/utest_constexpr_cmath.h @@ -1468,6 +1468,117 @@ bool test_fminf_rt() { return allCorrect; } +template +constexpr bool test_log_ct() { + static_assert(cxp::log::f(static_cast(1.0)) == static_cast(0.0), "log(1) should be zero"); + static_assert(cxp::log::f(static_cast(0.0)) == -std::numeric_limits::infinity(), + "log(0) should be negative infinity"); + static_assert(cxp::isinf::f(cxp::log::f(static_cast(0.0))), "log(0) should be infinite"); + static_assert(cxp::isnan::f(cxp::log::f(static_cast(-1.0))), "log(negative) should be NaN"); + static_assert(cxp::isinf::f(cxp::log::f(std::numeric_limits::infinity())), "log(infinity) should be infinity"); + static_assert(cxp::isnan::f(cxp::log::f(std::numeric_limits::quiet_NaN())), "log(NaN) should be NaN"); + + static_assert(cxp::log::f(static_cast(2.718281828459045)) > static_cast(0.9999) && + cxp::log::f(static_cast(2.718281828459045)) < static_cast(1.0001), + "log(e) should be approximately one"); + static_assert(cxp::log::f(static_cast(0.5)) > static_cast(-0.6932) && + cxp::log::f(static_cast(0.5)) < static_cast(-0.6931), + "log(0.5) precision error"); + return true; +} + +constexpr bool test_pow_ct() { + static_assert(cxp::pow::f(2.0f, 3.0f) > 7.0f && cxp::pow::f(2.0f, 3.0f) < 9.0f, + "pow(float) integer exponent failed"); + static_assert(cxp::pow::f(2.0, 3.0) > 7.0 && cxp::pow::f(2.0, 3.0) < 9.0, + "pow(double) integer exponent failed"); + static_assert(cxp::pow::f(9.0f, 0.5f) > 2.999f && cxp::pow::f(9.0f, 0.5f) < 3.001f, + "pow(float) fractional exponent failed"); + static_assert(cxp::pow::f(9.0, 0.5) > 2.999 && cxp::pow::f(9.0, 0.5) < 3.001, + "pow(double) fractional exponent failed"); + static_assert(cxp::pow::f(-2.0f, 3.0f) > -9.0f && cxp::pow::f(-2.0f, 3.0f) < -7.0f, + "pow(float) negative odd exponent failed"); + static_assert(cxp::pow::f(-2.0, 4.0) > 14.0 && cxp::pow::f(-2.0, 4.0) < 18.0, + "pow(double) negative even exponent failed"); + static_assert(cxp::pow::f(0.0f, 2.0f) == 0.0f, "pow(float) zero positive exponent failed"); + static_assert(cxp::isinf::f(cxp::pow::f(0.0, -1.0)), "pow(double) zero negative exponent should be infinity"); + static_assert(cxp::isnan::f(cxp::pow::f(-2.0f, 0.5f)), "pow(float) invalid negative base should be NaN"); + static_assert(cxp::isnan::f(cxp::pow::f(std::numeric_limits::quiet_NaN(), 2.0)), + "pow(double) NaN input should be NaN"); + return true; +} + +template +bool test_log_rt() { + bool allCorrect{true}; + const T tolerance = std::is_same_v ? static_cast(1e-5) : static_cast(1e-12); + + auto check = [&](T value) { + const T expected = std::log(value); + const T actual = cxp::log::f(value); + if (std::isnan(expected) ? !std::isnan(actual) + : std::isinf(expected) ? (!std::isinf(actual) || std::signbit(expected) != std::signbit(actual)) + : std::abs(actual - expected) > tolerance * std::max(1, std::abs(expected))) { + std::cout << "Failed: cxp::log::f(" << value << ") expected " << expected << " but got " << actual + << std::endl; + allCorrect = false; + } + }; + + check(static_cast(0.125)); + check(static_cast(0.5)); + check(static_cast(1.0)); + check(static_cast(2.0)); + check(static_cast(10.0)); + check(static_cast(-1.0)); + check(std::numeric_limits::infinity()); + + if constexpr (std::is_same_v) { + const float actual = cxp::logf::f(2.0f); + if (std::abs(actual - std::log(2.0f)) > 1e-5f) { + std::cout << "Failed: cxp::logf::f(2.0f) should match std::log(2.0f)" << std::endl; + allCorrect = false; + } + } + + return allCorrect; +} + +template +bool test_pow_rt() { + bool allCorrect{true}; + const T tolerance = std::is_same_v ? static_cast(2e-5) : static_cast(1e-12); + + auto check = [&](T base, T exponent) { + const T expected = std::pow(base, exponent); + const T actual = cxp::pow::f(base, exponent); + if (std::isnan(expected) ? !std::isnan(actual) + : std::isinf(expected) ? (!std::isinf(actual) || std::signbit(expected) != std::signbit(actual)) + : std::abs(actual - expected) > tolerance * std::max(1, std::abs(expected))) { + std::cout << "Failed: cxp::pow::f(" << base << ", " << exponent << ") expected " << expected + << " but got " << actual << std::endl; + allCorrect = false; + } + }; + + check(static_cast(2.0), static_cast(3.0)); + check(static_cast(9.0), static_cast(0.5)); + check(static_cast(-2.0), static_cast(3.0)); + check(static_cast(-2.0), static_cast(4.0)); + check(static_cast(0.0), static_cast(2.0)); + check(static_cast(10.0), static_cast(-2.0)); + + if constexpr (std::is_same_v) { + const float actual = cxp::powf::f(2.0f, 3.0f); + if (actual != std::pow(2.0f, 3.0f)) { + std::cout << "Failed: cxp::powf::f(2.0f, 3.0f) should match std::pow(2.0f, 3.0f)" << std::endl; + allCorrect = false; + } + } + + return allCorrect; +} + // Runtime tests to complement compile-time tests bool runtime_tests() { bool allCorrect{true}; @@ -1550,6 +1661,12 @@ bool runtime_tests() { allCorrect &= test_fmaxf_rt(); allCorrect &= test_fminf_rt(); + // Test log and pow with runtime values + allCorrect &= test_log_rt(); + allCorrect &= test_log_rt(); + allCorrect &= test_pow_rt(); + allCorrect &= test_pow_rt(); + // Test signbit with runtime values, where the intrinsic path is taken allCorrect &= test_signbit_rt(); allCorrect &= test_signbit_rt(); @@ -1615,6 +1732,9 @@ int launch() { static_assert(test_fminf_ct(), "fminf compile-time tests failed"); static_assert(test_fmax_fmin_double_ct(), "fmax/fmin double compile-time tests failed"); + static_assert(test_log_ct(), "log compile-time tests failed for float"); + static_assert(test_log_ct(), "log compile-time tests failed for double"); + static_assert(test_pow_ct(), "pow compile-time tests failed"); // Runtime tests if (!runtime_tests()) { diff --git a/utests/core/constexpr_libs/utest_constexpr_expf_exact.h b/utests/core/constexpr_libs/utest_constexpr_expf_exact.h index e8aa8f51..5cedf2fd 100644 --- a/utests/core/constexpr_libs/utest_constexpr_expf_exact.h +++ b/utests/core/constexpr_libs/utest_constexpr_expf_exact.h @@ -20,6 +20,7 @@ #include +#include #include #include #include @@ -83,6 +84,53 @@ int launch() { allCorrect = false; } + // Runtime dispatch check: outside constant evaluation cxp::expf must behave exactly + // like std::expf on a sweep of the useful input range, including boundaries. + { + constexpr float sweepStart = -105.0f; + constexpr float sweepEnd = 89.0f; + constexpr int sweepCount = 1001; + for (int i = 0; i < sweepCount; ++i) { + const float x = sweepStart + (sweepEnd - sweepStart) * static_cast(i) / (sweepCount - 1); + const uint expected = cxp::bit_cast(std::exp(x)); + const uint actual = cxp::bit_cast(cxp::expf::f(x)); + if (expected != actual) { + std::cout << "Runtime Fail: cxp::expf::f(" << x << ") expected bits 0x" << std::hex << expected + << ", got 0x" << actual << std::dec << std::endl; + allCorrect = false; + } + } + } + + // Compile-time sweep: values produced by the constexpr path are compared at runtime + // against std::expf. exp results are always non-negative, so the bit patterns are + // directly ordered and their difference is the distance in ulps (this also works + // for subnormals and for the 0/inf saturation points). + { + constexpr float sweepStart = -103.0f; + constexpr float sweepEnd = 88.7f; + constexpr int sweepCount = 257; + constexpr auto ctResults = []() { + std::array results{}; + for (int i = 0; i < sweepCount; ++i) { + const float x = sweepStart + (sweepEnd - sweepStart) * static_cast(i) / (sweepCount - 1); + results[i] = cxp::expf::f(x); + } + return results; + }(); + for (int i = 0; i < sweepCount; ++i) { + const float x = sweepStart + (sweepEnd - sweepStart) * static_cast(i) / (sweepCount - 1); + const uint expected = cxp::bit_cast(std::exp(x)); + const uint actual = cxp::bit_cast(ctResults[i]); + const uint ulpDiff = expected > actual ? expected - actual : actual - expected; + if (ulpDiff > 2u) { + std::cout << "Constexpr Fail: cxp::expf::f(" << x << ") expected bits 0x" << std::hex << expected + << ", got 0x" << actual << std::dec << " (" << ulpDiff << " ulps)" << std::endl; + allCorrect = false; + } + } + } + if (allCorrect) { std::cout << "All tests passed!" << std::endl; return 0;